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Abstract: We consider the problem of estimating an arbitrary smooth functional of k > 1 
distribution functions (d.f.s.) in terms of random samples from them. The natural estimate 
replaces the d.f.s by their empirical d.f.s. Its bias is generally ~ n~ , where n is the minimum 
sample size, with a pth order iterative estimate of bias ~ n~ p for any p. For p < 4, we 
give an explicit estimate in terms of the first 2p — 2 von Mises derivatives of the functional 
evaluated at the empirical d.f.s. These may be used to obtain unbiased estimates, where 
these exist and are of known form in terms of the sample sizes; our form for such unbiased 
estimates is much simpler than that obtained using polykays and tables of the symmetric 
functions. Examples include functions of a mean vector (such as the ratio of two means 
and the inverse of a mean), standard deviation, correlation, return times and exceedances. 
These pth order estimates require only ~ n calculations. This is in sharp contrast with 
computationally intensive bias reduction methods such as the pth. order bootstrap and 
jackknife, which require ~ n v calculations. 
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1 Introduction and summary 

Let T{F) be any smooth functional of one or more unknown distributions F based on 
random samples from them. Bias reduction of estimates of T(F), say T(F), has been a 
subject of considerable interest. Traditionally bias reduction has been based on well known 
resampling methods like bootstrapping and jackknifing in nonparametric settings, see Gray 
and Schucany (1972) and Efron (1982). However, these methods may not be effective in 
complex situations when the sampling distribution of the statistic changes too abruptly 
with the parameter, or when this distribution is very skewed and has heavy tails. Also the 
robustness properties of F may not be preserved for T(F) for all T(-). 

Recently, various analytical methods have been developed for bias reduction in paramet- 
ric settings. Withers (1987) developed methods for bias reduction based on Taylor series 
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expansions. Sen (1988) established asymptotic normality of y/n{T(F) — T(F)} as n — > oo 
under suitable regularity conditions. Cabrera and Fernholz (1999, 2004) denned a target 
estimator: for a given T and a parametric family of distributions it is denned by setting the 
expected value of the statistic equal to the observed value. Cabrera and Fernholz (1999, 
2004) established under suitable regularity conditions that the target estimator has smaller 
bias and mean squared error than the original estimator. See also Fernholz (2001). 

Suppose we have k > 1 independent samples of sizes n%, . . . ,n k from distribution func- 
tions (d.f.s) F = (Fi,...,F k ) on R Sl , . . . , R Sk . Let F = (Fi,...,F k ) denote their sample 
d.f.s and let n be the minimum sample size. The problem we consider in this paper is that 
of finding an estimate of low bias for an arbitrary smooth functional T(F). The natural 
estimate T(F) generally has bias ~ n , that is, 0(n -1 ) asn-> oo. 

This paper has already been cited as an unpublished technical report in Withers and 
Nadarajah (2008). The estimates proposed here have been compared to alternatives. We 
showed in particular that our estimates consistently outperform bootstrapping, jackknifing 
and those due to Sen (1988) and Cabrera and Fernholz (1999, 2004). We also provided 
computer programs in MAPLE for implementation of the proposed estimates. 

The emphasis of this paper is to describe how to find estimates of low bias for T(F). 
Because of the material in Withers and Nadarajah (2008), the emphasis here will not be on 
numerical illustrations or applications. 

For the reader's convenience, in Section 2, we repeat the definition of functional deriva- 
tives and rules for obtaining them given in Withers (1988). In Section 3, we have a formal 
asymptotic expansion of the form 

oo 

ET{F) = Y,n' r C r , (1.1) 

r=0 

where Co = T(F). The coefficient of n~ r in ET(F), C r (F,T) = C r may be written in 
terms of the (functional or von Mises) derivatives of T(F) of order < 2r, and is given there 
explicitly for r < 4. 

From (jl.ip if a functional T( n ) (F) can be expanded as 

oo 

T H =^2n~^(F) 

i=0 

then 

oo 

ET (n) (f) =^n" i C l (T), 

i=0 

where 

% 

C i {T)=Y i C i {F,T ir _ i ). 
i=o 

Defining Tj iteratively by To = T and 

i 

Ti(F) = -^2C j (F,T i _ j ) (1.2) 
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for i > 1 it follows that for p > 1 

p-i 

Tn P (F) = ^n" i T i (F) (1.3) 
i=0 

satisfies a formal expansion of the form 

ET np {F) = T(F) + O {n- p ) . 

So, T np (F) is a pi/i order estimate in the sense that it has bias 0(n _p ). This result was 
given for the case k = 1, p = 2 using a different approach in an unpublished technical report 
by Jaeckel (1972). 

Note that Ti{F) given by (jl.2p is the coefficient of n~ % in the expansion in powers of 
n^ 1 of the unbiased estimate (UE) of T(F), if an UE exists. 

Section 4 gives Ti{F) explicitly in terms of the first 2i derivatives of T(F) for i < 3. So, 
T n i{F) is an explicit estimate of bias 0(n -4 ). Note 4.1 shows how to obtain from fjl .3j) an 
estimate of bias 0(n~ p ) of the form S np (F), where 

p-i 

S np (F) = J2S l (F)/{(n-l)...(n-i)}. 

8=0 

This estimate is unbiased for one sample if T(F) is a polynomial in F (such as a moment 
or cumulant) of degree up to p. 

Section 5 gives examples and makes comparisons with the UEs of central moments and 
cumulants given by James (1958) and by Fisher (1929). Our method is demonstrated to 
be give much simpler results for UEs of products of moments than the polykay system of 
Wishart (1952) as expounded in Section 12.22 of Stuart and Ord (1987) using tables of the 
symmetric functions. 

Examples 5.1 to 5.3 estimate an arbitrary function of the vector n(F), the mean of one 
multivariate distribution. Example 5.2 specializes to T(F) = &' (i(F) /h' fi(F) , where a, b 
are given si-vectors, in particular for the ratio of means of a bivariate sample, 

T(F)= f x 1 (F)/ f x 2 (F). 

Examples 5.4 and 5.5 estimate an arbitrary function of the means of k univariate distri- 
butions; in particular it considers the case of two univariate samples (k = 2, s± = s 2 = 1) 
with 

r(F) = /i(F 1 )/ / u(F 2 ). 

Example 5.6 gives an explicit expression for the general derivative of the rth central moment 
(jl t . Together with the chain rule of Appendix A this enables one to obtain a pth order 
estimate of any smooth function of moments. In particular, we give fourth order estimates 
for any central moment and UEs for fj, r for r < 7. 

Examples 5.7 to 5.11 extend this to an arbitrary product of moments. An alternative 
matrix method for obtaining UEs of products of moments is given there. This involves 
obtaining simultaneously the UEs of all moment products of a given degree. Examples 5.12 
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to 5.15 give fourth order estimates of the standard deviation and functions of it. Example 
5.16 gives third order estimates of the ratio of the mean to the standard deviation. 

Examples 5.17 to 5.21 give applications to return times and exceedances. Examples 
5.22 and 5.23 illustrate how to obtain UEs for multivariate moments and cumulants from 
univariate analogs. Finally, Examples 5.24 and 5.25 give second order estimates for the 
correlation and its square. 

The method can also be used to estimate with reduced bias any cumulant of T(F). 
This is illustrated in Section 6 which gives a third order estimate for the covariance of any 
estimate of the form T(i ? ), where now T may be a vector. For example, by Example 5.1, 
if k = 1 and T(F) is any function of n(F) (such as (J*i(F) / ^{F)) if s± = 2), this estimate 
is a function of the mean and covariance of F only, whereas C% depends also on the third 
moment. 

Section 7 shows how to estimate the covariance of an estimate of bias. 

There are, of course, other pth order estimates of T(F), but they are all computationally 
intensive, requiring 0(n p ) calculations (except in special cases), whereas our method requires 
only 0(n) calculations for fixed p. The main examples are, firstly, the (p — l)th iterated 
bootstrap, 6 p -i of equation (1.35) of Hall (1992) in which (— 1) ,+1 should be inserted in the 
right hand side; and, secondly, the pth order jackknife of equation (4.17) of Schucany 
et al. (1971), a ratio of p x p determinants. To see that this requires 0(n p ) calculations 
note that t p of their equation (4.19) requires 0(n p ) calculations. 

The techniques given here can also be applied to quantify their biases. Note that if A 
and B are two pth order estimates of T{F) then A — B = O p (n~ p ). 

Appendix A gives a very useful chain rule for obtaining the derivatives of a function of a 
functional. Appendix B gives some results used to obtain {Tj} of fjl .3|) . Appendix C shows 
how to estimate the number of simulated samples needed to estimate the bias to within a 
given relative error. 

Tiit (1988) by an entirely different method obtained an expansion of the form (jl.ip for 

s 

m(v) = T(F) = Y[EX V % 
i=i 

where X ~ F, and so also for fi r (F). For these cases he constructs estimates of bias 0{n~ p ) 
given p > 1. He shows for T(F) = m(v) that the UE T noo (F) converges if E\X\ h < oo, 
where h = Y2t=i v i an d n — 1 > the number of partitions of h. His expression on page 12, 
Theorem 4, is incorrect. He gives 

var rh(v) = n~ l V + O {n 2 } , 

where 

s 

V = m (v) 2 [A — s 2 ) , and A = m2 Vi m 2 .. 

i=l 

Here, A should be 

s 

m Vi+V] m-^m-f. 
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For the case T(F) = ^ 3 his Table 2 illustrates through simulations for F = £7(0,1) and 
n = 5, 10 how the bias of T np (F) falls to zero as p increases. 

2 Functional partial derivatives and notation 

Let F s denote the space of d.f.s on R s . Let x, y, xi, . . . ,x r be points in R s , F G F s and 
T : F s — > R. In Withers (1983) the r-th order functional derivative of T(F) at (xi, . . . ,x r ) 

Jx lr .,x r = Tp (xi, . . . ,x r ) , 

was defined. It is characterized by the formal functional Taylor series expansion: for G in 

oo „ r r 

T(G) - T(F) 2> . . . , x,.) IJ d (G - F (x;)) /r!, (2.1) 

r=l 7 j=l 

where J r denote r integral signs, and the constraints T Xl ... Xr is symmetric in its r arguments, 
and 

J T^.^dF (xi) = 0. 

These imply F(x.j) in (|2.ip can be replaced by zero. In particular, it was shown that, for 
< e < 1, 

T« = 9T(F + e((y ai -F))/Se 

at e = 0, where <5 X is the d.f putting mass 1 at x, that is S x (y) = I(x < y) = 1 if x < y and 
otherwise. For example, T(F) = F(y) has first derivative T x = Tp{x) = 5 x {y) — F(y) = 
F{y) x , say. 

Also, T Xl ... Xr = if T(F) is a 'polynomial in F : of degree less than r (for example, a 
moment or cumulant of F of order less than r), so that the Taylor series in f)2. 1 1) consists 
of only r — 1 terms. Note that T(F) is a polynomial in F of degree m if for any G in J- s , 
T(F + e(G — F)) is a polynomial in e of degree m. 

Suppose now that F = (Fi, . . . , F^) consists of k distributions on R Sl ,...,R Sk and 
that T(F) is a real functional of F. Then the functional partial derivative of T(F) at 

( Ql Qr | is defined by 

V xi • • • x r J 

T^ 7=T F ( ai "' ° r V 

1 r v x i • • • X W 

where Xj in iZ 5 "* and Oj in {1, 2, . . . , A;}, and is obtained by treating the lower order functional 
partial derivatives and T(F) as functionals of F a alone for a = a\,...,a r . For example, 
^xi-'-xp ^ s ^ ne ordinary functional derivative of S(F a ) = T(F) at (xj • • • x r ), and ^xi-.^yi-ys 
is the ordinary functional derivative of S(Fb) = T®"'.? x at (yi ■ • • y s ). 

Just as d 2 f (x , y) / dxdy = d 2 f(x,y)/dydx under mild conditions, swapping columns of 
T ^ °^ (for example, ^ and x ^) will not alter its value. So, r x ^'.'. a x ^y 1 ... ys is also the 
ordinary functional derivative of S(F a ) = Ty^ yr at (xi • • • x s ). 
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The partial derivatives may also be characterized by the formal functional Taylor series 
expansion: for G = (G\ ,...,(?&) £ F Sl x • • • x F Sk , 

00 pr / \ r 

T(G) - T(F) « £ / T ? x! x II d ( x i) " ^ M M (2-2) 

r=l ^ ^ r ' j=l 

with summation of the repeated subscripts a\ ■ ■ ■ a r over their range 1 ■ ■ ■ p implicit, together 
with the constraints 

^xi - xr * s n °t altered by swapping columns, 

and 



/ T%::%dF ai ( Xl ) = 0. 



These imply F aj (xj) in (|2.2p can be replaced by zero. The partial derivatives may be 
calculated using 

T» = 5 X for S(F a )=T(F), (2.3) 

and 

ctl • • • a r +i 



T 

xi • • • x r+ i 



Ol ■ ■ ■ a r \ a r+ i ST~^ c ^ /a r --a r+ i \ /r) /|\ 

Jil ■ ■ ■ -*-r / — ' 
1=1 



where <5jj = 1 or for i = j or i / j and ()j means 'drop the i-th column'. For example, 
+ 5 ab T^. The proof of $£$-$23$ is as for equation (2.6) of Withers (1983). 



3 Expansions for bias 

Perhaps the easiest method to obtain expressions for the bias coefficients {C r } of (jl.ip and 
the bias reduction coefficients {Ti(F)} of (II. 3D is from their parametric analogs, given in 
equation (A.l) and Appendix D (for i < 3) of Withers (1987). The method is to identify 
(0,0,t,^2) with (F,F,T, J), where the integral is with respect to the appropriate distri- 
bution F{. This method was used in Withers (1988) to derive non-parametric confidence 
intervals of level 1 — a + 0(n~i/ 2 ) from their parametric analogs. It is convenient to set 

T(aV ...) /•••/'//( % •••) dF a (x)dF b (y) • • • , (3.1) 

where x l denotes a string of % x's (not a product) and similarly, for a\ In the notation of 
Withers (1988) this is [1*2* • ••]„&.•- Setting 

A a = n/n a with n = min rij, (3.2) 

the above approach yields 

Cx = \2\/2, C 2 = |3|/6 + |2 2 |/8, (3.3) 
C 3 = |4|/24 + |23|/12 + |2 3 |/48, (3.4) 
C 4 = |5|/120 + |24|/48 + |3 2 | /72 + j 2 2 3 j /48 + |2 4 | /384, (3.5) 
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where 

|2| = ^A a T (a 2 ), 

|3| = E A ' T (« 3 )' 

|2 2 | =^A ai A a2 T(a 2 a 2 ), 

|4i = E A '{ r (« 4 )- 3T (« 2 « 2 )}' 

|23|=^A a A 2 T (a 2 6 3 ), 

|2 3 | =^A ai A a2 A a3 T(a 2 a 2 a 2 ), 

|5i = E A '{ T (« 5 )- 10T ( aV )}' 

|24|=EAaAf{T(a 2 6 4 )-3T(aW)}, 

|3 2 |=£« 2 r (^i)> 

|2 2 3| =E A a 1 Aa 2 A^T (a 2 a 2 . 6 3 ) , 
|2 4 | = A ai A a2 A a3 A a4 T (a 2 a 2 a 2 a|) . 
For example, if k = 1 (one sample) then 

d=T (l 2 ) /2, C 2 = T (l 3 ) /6 + T (l 2 l 2 ) /8, . . . . (3.6) 

More generally, 

\A i \=E<- 1 ---<- 1 \ Ai l 1 -a^ 

\**\ = E A^ 1 • • • A^Ag- 1 • • • Ag" 1 \A*B>\ ai ^ bj (3.7) 

with each a±, . . . ,bj summed over 1, . . . , k, 

\A^\ ai ^ bi bj =T(af... afbf • • • bf) if A and B = 2 or 3, 

|4| a = r(a 4 )-3T(aV), 
\5\ a = T (a 5 ) - WT (a 2 a 3 ), 
|24| a6 = T(a 2 6 4 ) - 3T (a 2 6 2 6 2 ) . 



For example, 



and 



a\a,2 



\A 2 \ = T(a A a A ) if A = 2 or 3 

I 1 0.10,2 * ' 

= JJt f ($ f A yF ai (x)dF a2 (y), 

so for the one sample case (k = 1), 

\A'\ =T(l A ---l A ) if ,4 = 2 or 3, 

|A*J3 J "| = T (l A • • • • • • I s ) if A and 5 = 2 or 3, 

|4| = T (l 4 ) - 3T (l 2 ! 2 ) , |5| = T (l 5 ) - 10T (l 2 ! 3 ) , |24| = T (l 2 ! 4 ) - 3T (l 2 l 2 l 2 ) . 
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Note 3.1 The general term C r is given by equation (A.l) of Withers (1987), \3.2) , {3. 7\ ), 
and 

\ij ■ ■ ■ lab- = J d l K a (xi • • • J dK' b (yi • • • yj) • • • 7> f ^ _ ^ ^ • • • V 

where J 1 d l K' a (xi • • • Xj) is £/ie Lebesgue-Stieltjes integral, 

xi A X2 A • • • = min (xi, X2, . . .) taken componentwise , 
/l2- = -Fa (xi A x 2 A • • • ) , 

K a (xiX2 • • • ) = k (Yi, Y 2 , . . .) , t/ie joint cumulant at Yj = I (X a < Xj) , 
«4 (xix 2 • • • ) = K a (xix 2 • • • ) expressed as a function of {fij—} at /j = 0, 

and I is the indicator function and X a ~ F a . For example, using an obvious summation 
notation 

K a (xix 2 ) = fu - /1/2, 
3 

K a (X1X2X3) = /123 - ^ /12/3 + 2/1/2/3, 

4 3 
K a (xi • • • X 4 ) = /1...4 - ^ /123/4 - ^ /l2/34, 

3 

«1 (X1X 2 ) = fi2, K a (xiX 2 X 3 ) = /123, K' a (Xl, . . . ,X 4 ) = /1...4 - ^ fizfai- 

Note 3.2 As a check if k = 1, (Ci,C 2 ) = (an, 012) on page 580 of Withers (1983). 

4 Estimates of bias 0(n 4 ) 

Here, we give expressions for {Tj,i < 3} of (jl.2p and for {<Si,i < 3} of Note 4.1. Estimates 
of bias 0(n~ 4 ) are then given by T nA [F) of (JUSJ) and 5 n4 (F) of (03]), M . 

From their parametric analogs in Appendix D of Withers (1987), we obtain (see Ap- 
pendix B) in the notation of (|3. T[) 

Ti(F) = -|2|/2, T 2 (F) = |3|/3+|2 2 |/8-^A2T(a 2 )/2, (4.1) 

and 

T 3 (F) = -^A^(a 2 )/2 + ^A3r(a 3 )-^A^(a 4 )/4 + ^A3r(aV)/2 

+ 2 A a A 2 T (a 2 b 2 ) /4 - J] AaA fe 2 T (a 2 6 3 ) /6 - J] A a A b A c T (aW) /48. 

For the one sample case (A; = 1), these reduce to 

T!(F) = -T(l 2 )/2, (4.2) 

T 2 (F) = T (l 3 ) /3 + T (1 2 1 2 ) /8 — T (l 2 ) /2, (4.3) 
T 3 (F) = — T (l 2 ) /2 + T (l 3 ) - T (l 4 ) /4 + 3T (l 2 l 2 ) /4 — T (l 2 l 3 ) /6 

-T(l 2 l 2 l 2 ) /48. (4.4) 
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Note 4.1 Let {Ni(n),i > 0} be given functions satisfying Ni(n)/n 1 — > 1. Then U.3\) may 
be rewritten as S np (F) + 0(n~ p ), where 

p-i 

S np (F) = ^2N i (n)S i (F). (4.5) 

So, S np (F) is a pth order estimate ofT(F). 

Suppose now that it is known that there exists an UE and that it has the form S np {F). 
Then this gives a method of obtaining it. For example, if k = 1 and T(F) is a polynomial 
of degree p in F (for example, a product of moments or cumulants of total degree p ), then 
the UE of T(F) has the form g3|) with 

N i (n) = l/(n-l) i , (4.6) 

where (r)j = r!/(r — i)! = r(r — 1) • • • (r — i + 1). In this case, {Si} are given in terms of 
{Ti} by equation (2.11.2) of Withers (1987): 

So = T, Si = Ti, 52 = T 2 — T\, S3 = T3 — 3T2 + 2Ti, .... 

If k = 1 and we choose Ni(n) as in |^.6[ ), then Si is generally a simpler expression than Ti: 

S (F)=T(F), S 1 (F) = -T (l 2 ) /2, 

S 2 (F) =T(l 3 )/3 + T(l 2 l 2 )/8, (4.7) 
S 3 (F) = -T (l 4 ) /4 + 3T (1 2 1 2 ) /8 - T (l 2 l 3 ) /6 — T (l 2 l 2 l 2 ) /48. (4.8) 

//Ml, 

5 (F) = T(F), 5i(F) = 0/ J£f , 

S 2 (F) = T 2 (F) - T\{F) = |3|/3 + |2 2 | /8 + ]T (A a - A 2 ) T (a 2 ) /2, 

and so on. 



Note 4.2 Forp > 1, sei e HtP (T,F) = T np (F) of i TOI) and Zei {^(F)} 6e smooth. T/ien a 
pth order estimate of 

00 

U n (F) = Y j n- i U i {F) 

i=0 

is 

p-i 

U tn) P (F) = E ™~ (^) • ( 4 - 9 ) 

8=0 

Zei « r (X) denote any rth order cumulant o/X, any qxl random vector. Then n 1_r K r (T(F)) 
can be expanded in the form |^.g[ ); a method of obtaining {Ui} is illustrated in Section 6 for 
the case r = 2. 
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Note 4.3 ET(F) may be infinite or may not exist. For example, this is the case if k = 
s = 1, T(F) = , I > 1 and F has positive density at zero, or F(x) approaches zero 

too slowly as x — )■ 0. So, Quenouille (1956, page 356) is wrong in giving X finite bias for 
X ~ N (2,1). In such cases, our method may be salvaged provided we know an upper bound 
for \T{F)\, say \T[F)\ < u < oo. By large deviation theory P(\T(F)\ > u) = 0(exp(— n\)), 
where A > 0. Typically, T np (F) is a pth order estimate of'T(F), where 

f np (F) = \ T ^ (4.10) 
[ c, otherwise 

and c is any finite constant, for example, u. 

The estimates J^.5| ) and \4-9\ ) can be adapted similarly, to give S np (F) and U* p (F) say. 
Similarly, if U^(F) is the formal expansion of n r ~ 1 K r (T np (F)) then 

U* q (F)I ^ T(F) < u j is a qth order estimate of n r ~ l n r \ T np (F) 

even if n r (T(F)) is not finite. For example, the variances in equations (10.17)-(10.20) of 
Kendall and Stuart (1911) are infinite if the density at zero is positive. 

Note 4.4 An alternative estimate of bias 0(n~ p ) is T^ p {F) = T nq {F), where q < p is the 
maximum integer such that {n~ l Ti(F),0 < i < q} decreases in absolute value. This may be 
useful if T np (F) diverges. Note that S^ p (F) and T^ (F) may be defined analogously from 
ggP and ffTlfy . 



5 Examples 

Example 5.1 Suppose k = 1, X ~ F on R s and T{F) = g(fi), where fi = n{F) = EX. has 
dimension s\ = s and g is a function with finite derivatives at fi. By the chain rule iA.6\ ) 
or (A.l\ ) of Appendix A, 

Tf ( X l) • • • , X r ) = gj 1 ---j r fJ-j 1 xi ' ' ' HjrlLrl 

where 

^'x = Mjf(x) = Xj - fij, g ■ ■ ■ = g ■ ■ ■ (n) 

are the partial derivatives of g(fi) with respect to \i, and summation of the repeated indices 
ji, . . . ,j r over their range 1, . . . , s is implicit. So, 

T (l n l i2 • • • ) = ..,../.,. ../,,,...// [ji • ■ -j h ] n [A* • • • k i2 ] 

where 



b'l •••3a] = J (xh - Mil) ■ • • ( x ja ~ Mia) dF(x), (5.1) 
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the joint central moment. So, 

s 

T (l 2 ) = 9ijfJ-[ij] = ^2gufi[n] + 2 9ijV[ij], 
i=l i<*<i<s 

T (! 3 ) = 9ijk^[ijk], 

T (! 4 ) = 9ijkiv[ijkl], 

T (l 2 l 2 ) = gj lj2kl k 2 ^ [jik] M [hfa] , 

T(1 2 1 3 ) = gijki m v[ij]n[klm], 

T (1 2 1 2 1 2 ) = g ijMmn n[ij}n[kl]n[mn\. 

So, by (PP-BP 

T 1 (F) = -C 1 = -g ij n[ij}/2, 

T 2 (F) = -g ij n[ij}/2 + g ijk n[ijk]/3 + g ijk ifi[ij]fi[kl]/8, 
T 3 (F) = -g ij fi[ij}/2 + gijkfJ>[ijk] - g ijk i {^[ijkl] - 3fi[ij]fi[kl]} /4 
-9ijklmlJ>[ij]tJ>[klm]/6 - g ijk i mn fi[ij}ii[kl}ii[mn}/48. 

A pth order estimate ofT(F) is now given in terms of these by T np (F) of tl.3\) . 

Example 5.2 Consider Example 5.1 with g{n) = a'fi/(3'fi = N/D, say, where a, j3 are 
given s-vectors. Its ith order partial derivative with respect to n is 

i 

g h ... jt = {-l) l -\i - l)\D-* Sjifa ■■■Pov (5-2) 

where 

Si = Oi- PiT(F) (5.3) 

and 

in 

^ ] fil—i m fil—im fi2—imil + " " " + fimH—im-1' 

So, 

t (r) = ■ ■ ■ fan [n ■ ..ji] , 

t (i 2 i 2 ) = v.n % x fafafav Uih] » [kk] , 

T (1 2 1 3 ) = AID' 5 {25 n /fa + 36 js /fa}fa • • • fan [hj 2 ] n lhhj 5 } , 

t (i 2 i 2 i 2 ) = -6\D- e s h fa ■ ■ ■ fan [j x j 2 ] n [kk] /* Ikk] ■ 

In particular, for g(fi) = nil 1 V>i (the ratio of means for one bivariate sample), 
T (l<) = (-ly-Hlrf [12 4 - 1 ] - T(F)n [?] } , 

t (i 2 i 2 ) = -v.nt [2 2 ] - T{F)n [2 2 ] 2 } , 

T (1 2 1 3 ) = 4!/zJ 6 {2n[l2]n [2 3 ] + 3/i [2 2 ] n [12 2 ] - 5T(F)n [2 2 ] n [2 3 ] } , 

t (i 2 i 2 i 2 ) = -einf {n[i2} - T(F)n [2 2 ] } n [2 2 ] 2 , 
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so 



S 1 (F) = Tx(F) = -d = ^ 2 {^[12] - T(F)fi [2 2 ] } , 

T 2 (F) = 2^ 3 [12 2 ] - T(F)fi [2 3 ] } — {l + 3/^ V [2 2 ] } , 

S2(F) is the same as T 2 (F) with '1 +' deleted, 

T 3 (F) = ^ 2 { M [12] - T(F)» [2 2 ] } {l - 18^ V [2 2 ] - 8/^ V [2 3 ] + [2 2 ] '} 

+6/^ 3 {/* [12 2 ] - T(F)fi [2 3 ] } {1 - 2/zJ 2 /z [2 2 ] } 
+6^ 4 { M [12 3 ] - T(F)^ [2 4 ]} 

cmo? 

5 3 (F) = iq 2 { M [12] - T(F)fi [2 2 ] } {-9^V [2 2 ] - [2 3 ] + 15/zJ V [2 2 f} 

-12^ 5 {/. [12 2 ] - T(F)n [2 3 ] } M [2 2 ] + 6^ 2 " 4 [l2 3 ] - T(F)» [2 4 ] } . 

[Table 5.1 about here.] 

Example 5.3 Consider Example 5.1 with g(fi) = (a'fi) p = N p , say, where a is a given 
s-vector. The ith order partial derivative of g(n) with respect to fi is 



Set 



Then 



9h-n = {PW P ''>.,: 

<*(*) =N~ % a jl ---a ji n\j 1 ---j i ]. 

T{l i ) = (p) t NPa (l) , 
T (1 2 1 2 ) = (p) 4 iV p a 2 2) , 

T(1 2 1 3 ) = (^5^0(2)0(3), 

T(1 2 1 2 1 2 ) = (p) 6 NPaf 2) , 
T 1 (F) = -C 1 = -(p) 2 N*a {2) /2, 

T 2 (F) = ^ {-(p) 2 a (2) /2 + (p) 3 a (3 )/3 + (p) 4 a 2 2) /s} , 
T 3 (F) = NP | - (p) 2 a (2) /2 + (p) 3 a (3) - (p) 4 [a (4) - 3a 2 2) ] /4 

3 



-(p) 5 a (2) a (3) /6 - (p) 6 a (2) /48 

In particular, for a univariate sample (s = 1) with central moments {n r } and g(fj,) = fi p , 

S 1 (F) = T 1 (F) = -(p) 2 ^- 2 f i 2 /2, 
T 2 {F) = -(p) 2 ^- 2 M2 /2 + S 2 (F), 
S 2 {F) = (ph^-^z/S + ip^-^j/S, 

Ts(F) = -( p ) 2 ^-V 2 /2 + (p) 3 ^-V3-(p)4/i p - 4 (^-3^)/4 
-(p) 5 /i p ~W 2 /6 - (p) 6 ^~Vi/48 
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and 

S 3 (F) = -Cp) 4 /i p - 4 (2 M4 - 3^1) /8 - (p) 5/ u p - W2/6 - (p) 6 ^- 6 /ii/48. 

In particular, for p a positive integer, by Note 4-1, an UE for [i v is 

p-i 

(F)/(n-l) f , 

where So(F) = [i v , and 

forp = 2:S x {F) = - l P 2 , 

forp = 3: S X (F) = -3^ 2 , S 2 (F) = 2 M3 , 

/or p = 4 : 5i(F) = -6/i 2 /i 2 , S 2 (F) = 8/4/3 + 3/i 2! 5 3 (F) = -6/x 4 + 9/4 

These results may be checked by by solving the system of equations given by Wishart (1952, 
page 5). Forp = 4 the system has seven equations. Alternatively, one may follow the method 
of Section 12.22 of Stuart and Ord (1987) using their tables of the symmetric functions. 
For example, after some labor one obtains for p = 4 the UE T n (F), where 

(n-l) 3 T n (F) = (N 3 - 8n 2 + 23n - 30) m 4 - n (n 2 - 7n + 4) m 3 m 1 
—n (n 2 — 6n + 6) m\ + n 2 (n — 9)m 2 mf + n 3 mf, 

where mi = EX 1 . Clearly, our method gives a much simpler form. 
For p = —1, that is T(F) = jjL~ l , the above gives 

p-i 

S np (F) = Y,S l (F)/(n-l) u 

where 

S (F) = fi- 1 , S 1 (F) = -n- 3 » 2 , 
S 2 (F) = -2/i- 4 /X3 + 3/i- 5 /x|, 

S 3 (F) = -3^~ 5 (2/z 4 - 3/4) + 20 A r f W 2 - 15aT 7 A& 
so setting j r = fi r fj,~ r , Sj = Si(F)/T(F) is given by 

si = -72, 

S 2 = -273 + 37|, 

s 3 = -3 ( 2 74 - 3 7 |) + 2O7372 - 15 72 3 . 

Some simulations estimating the bias of S n i(F) of fl^.5[ ), and Note 4-3 with c = 1/u = 

n/W for 1 < i < 4, for , are given in Table 5.1. For n = 100 and p = 2, the estimates 
are poor: see Appendix C. 

Example 5.1 estimated a smooth function of the mean of one multivariate distribution. 
We now estimate a smooth function of the means of k univariate distributions. 
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Example 5.4 Suppose we have k univariate samples (that is s% = ■ ■ ■ = s& = 1) with 
T(F) = g(p), where now /x = (fj,(F±), . . . , /u(Ffc)). That is, T(F) is a function of the means 
of k univariate samples. Then 

f ax ■ ■ ■ a r \ _ 

J-F I I — 9a\ - a r l^a\X\ ' ' ' l^a r x r i 



X\ ■ ■ ■ X 

where g ■ ■ ■ is the partial derivative with respect to fj, and 

Hax = HF a (x) = X - H (F a ) =X — H a . 

So, 

T (tfW ■■■) = <l l ri„...li; (iji ! h\ 

where 



Mi [a] = Mi ( F a) = J (x - Ma) 4 dF a (x), 
the ith central moment of F a . So, for X a of 113. fy) . 
Ci = A a g aa /x 2 [a]/2, 

a 

a ab 
°3 = Yl X a9aaaa {fJ>M ~ ^[a] 2 } /24 

+ Y ^a>%gaabbbH2[a}nz[b)i '12 + ^ A A 6 A c # ao66cc M[a]M2 [&]M2 [c]/48, 
7\(F) = -C l5 

32(F) = ^A^ aaa /i 3 [a]/3 + ^A a A b 5 aafefe /ii[a]/X2[6]/8-^A^ aa ^ 2 [a]/2, 
7s(F) = -^A^ aa/ i 2 [a]/2 + ^A^ aaa/ i 3 [a] -^A^ aaaa { M4 [a]/4 + /i 2 [a] 2 /2} 

+ ^2 AaA 6 5f aa6fc ^ 2 [a]/x 2 [6]/4 - ^ A a A^ aaW , fe /i 2 [a]// 3 [&]/6 

- A a A A c g aabbcc ^ 2 [a]^ 2 [6]^ 2 [c]/48. 

Example 5.5 Consider Example 5.4 with g(fj,) = a'n/f3'n = N/D, say, where a. and (3 
are given k-vectors. Set 

la = a a //3 a - T(F), 

a 

Biki = {Aiki} at 7 a = 1, 
Ak = A ki, 
Bk = Bofci- 
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Then, by I5Hj> . 



Ci = -A 2 , 

C 2 = A 3 -6A 2 B 2 , 

C 3 = -A 4 + 3A) 22 + 6A 2 B 3 + 9A 3 B 2 - 15A 2 Sf , 
Ti(F)=A 2 , 

r 2 (F) = 24j - 3A 2 £ 2 + A m , 

T 3 (F) = A 22 i - 9A 131 - 3A 3 + 6A 4 - 12A) 22 - 3A 121 B 2 - 3A 2 B 121 - 8A 2 B 3 
-Y1A 3 B 2 + 15A 2 Bl 

In particular, for g{fi) = f-ii/f-12 (the ratio of means for two univariate samples), setting 
Vk = H 2 k Hkfe], we obtain 

C\ = \ 2 u 2 m/^ 2 , 

C 2 = \ 2 2 (-f 3 + 6z/f) flx/fi 2 , 

C3 = A2 (va - 3i/| - 15^2^3 + 15f|) (J-l/fJ-2, 

Ti(F) = -\ 2 v 2 fi\/ fi 2 , 

T 2 (F) = X 2 2 (-21/3 - ^2 + 3^ 2 2 ) 

r 3 (^) = A2 (-61/4 - 6^ 3 - i/ 2 - 15^| + 20^2 + 18^ 2 ) Ml/M2- 
This may also be derived from A5.2\) . 

Central moments and functions of them may be viewed as functions of noncentral mo- 
ments and so dealt with using Examples 5.1 and 5.4. However, it is much more convenient 
to deal with them directly in terms of the derivatives of the central moments. We now give 
these. 



Example 5.6 One univariate sample ( that is k = s\ = 1) with T{F) = fi r (F) = fir> 
the rth central moment of X ~ F. Let [i = fi(F) denote the mean of F. Recall that 
(r)i = r!/(r — £)! and set hi = fi Xi = x% — fi. The general derivative of fi r (F) is 



v 



T Xl ... Xp = fi rF (xi • • • x p ) = (-l) p < (r) p fi r - p - (r) p _i ^ [h J p - /ir-p+ify *J > hj.(5A) 

{ i=l J j=l 

For example, 

T x = -rfi r -\fi x + ff x - fi r , 

2 

T xy = (r) 2 ft r - 2 fi x (i y - r ^ (/i^r 1 - Hr-i) fiy, 

xy 

3 

T xyz = -(r)3fl r -3lixflyflz + (r) 2 (^" 2 - /V-2) MyMz- 

xyz 

These basic building blocks are written out more explicitly up to r = 6 in Appendix D. 



15 



Setting q = ii + i% -\ , this gives 



")q^r—q /^i 



7=1 jV/ 
0, ifq>r, 

(-lrV-i)!]!^' ^ = r - 



For example, 



/i r (l 2 ) = (r) 2 /i r -2M2 - 2r/i r , (5.6) 

/x r (l 3 ) = -(r) 3 ^ r _3/i 3 + 3(r) 2 (/U r - ,u r _ 2 /i 2 ) , (5.7) 

/i r (l 4 ) = (r) 4 ^ r _ 4 /x 4 - 4(r) 3 (/i r - , (5.8) 

/i r (l 2 l 2 ) = (r) 4 ^ r _ 4 /X2 - 4(r) 3/ u r _ 2 /i 2 , (5.9) 

/i r (l 2 l 3 ) = -(r) 5 /i r _ 5/ u 3/ u 2 + (r) 4 (2/i r _ 3 /i 3 + 3/j, r _ 2 /"2 - 3^ r _ 4/ u 2 ) , (5.10) 

^ r (l 2 l 2 l 2 ) = (r)6Mr-6A*2 ~ 6(r) 5 ^ r _ 4/ U2- (5.11) 



Substituting into the expressions of /l3.3\) - {375\) for the coefficient Ci of n 1 in the expansion 
of E\x r {F) gives 

Ti(F) = -d = r/i r - (r) 2/ u r _ 2 /i 2 /2, (5.12) 
C2 = (r) 2 ^ r /2 - (r) 2 (r - l)/i r _ 2 ^ 2 /2 - (r) 3 ^ r _ 3/ u 3 /6 + (r) 4 ^ r _ 4 ^/8, (5.13) 
C3 = -(r)3Vr/6 + Mr-2M2(r) 3 (r - l)/4+ (r) 3 (r - 2)^ r _ 3 ^ 3 /6 

+ (r) 4 /i r _ 4 (/i 4 - 3(r - l)/u 2 ) /24 - (r) 5 ^ r _ 5/ u 3 /i 2 /12 

+ (r) 6 /i r _ 6 ^/48, 
C 4 = (r) 4 ^ r /24 - (r) 4 (r - 7)^ r _ 2/ u 2 /12 - (r) 6 ^ r _ 3 /i 3 /2 

+// r _ 4 {-(r) 4 (r - 3) ^ 4 /24 + (r) 4 (r 2 - 3r - 8)^ 2 /16} 

+Vr-5 {-(r) 5 ^ 5 /120 + (r) 6 (r - 2)^^/12} 

+(r) 6 /i r „ 6 (/i 4 ^ 2 /48 + ^|/72 - r^/48) - (r) 7 /v-7/^l/48 

+(r) 8 /i r _ 8 //|/384. 



Substituting into the expressions of ft4.3\ )- fi4~4\ ) for the coefficient Ti(F) of n 1 in the expan- 
sion for the UE of /i r (F) gives 

T 2 (F) = r 2 fi r - (r 3 - r) /i r __ 2 ^ 2 /2 - (r) 3 /i r _ 3 /^ 3 /3 + (r) 4 /i r _ 4 ^/8, 

and 

T 3 (F) = rV - (r 4 - r) ^r-2^/2 - (r) 3 (r + 3)/i r _ 3 /i 3 /3 

+(r) 4 )U r _ 4 {-2/i 4 + (r + 6)^2} /8 + (r) 5 /i r _5/i 3 /u 2 /6 - (r) 6 /i r _ 6 ^/48. 

Similarly, from J^. 7| ) and J^.ff[ ), 

5 2 (F) = (r) 2 /i r - r 2 (r - l)/i r _ 2 /i 2 /2 - (r) 3 ^ r _ 3) u 3 /3 + (r) 4 /i r _ 4 //|/8 
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and 

S 3 (F) = (r) 3/ u r - r(r) 3 ^ r _ 2 /i2/2 - r(r) 3 ^ r _ 3/ u 3 /3 

-(r) 4/ u r _ 4 /x 4 /4 + (4r - 9)(r) 4 /i r _ 4 ^2/8 + (r) 5 Mr-5M3/42/6 ~ (r) 6 fi r _ 6 (4/48. 

Now from James (1958, page 6) the UE for fj, r has the form 

l r = ^a ir (V W[](l-i/n) (5.14) 
I i=0 J i=l 

/or r = 2s or 2s + 1, which can be recovered from {Tt,i < s} as in Note 4-1- So, the above 
{Ti,i < 3} provide UEs for jj, r for r < 7. These were given for r < 6 on James (1958, page 
6) and agree with our results. 

For example, for T(l 2 ) = -2fx 2 , so S^F) = 3fi 3 and T(l 3 ) = 12/j, 3 , T(1 2 1 2 ) = 0, 
so 82(F) = 4/u 3 and so the UE o//u 3 is 

Ms (f) {1 + 3/(n - 1) + 4/(n - 1) 2 } = ^ (f) { (l - n" 1 ) (l - 2rT X ) j" 1 . 

For r = 7, we obtain in this way {a^ = ovK-F)} of ft5.14\ ) as 

a 7 = M7, «17 = -7 (2/x 7 + 3^ 5 /x 2 ) , 

a 27 = 7 (11^7 + 39/i 5 ^ 2 - 10^4^3 + 15^3/^1) > 

037 = -7 (28/z 7 + 192/i 5 ^ 2 - 80^ 4/ u 3 + 6O//3//2) . 

Example 5.7 One univariate sample (that is k = s\ = 1) with T(F) = Y\.^=2^/ f or {Pj} 
arbitrary and {/ij} as in Example 5.6. Set Si(fi) = ^ and g(S) = n^f*- The ordinary 
partial derivatives of g(S) are 

9i = PiVi l T{F), gij = pi (pj - 8ij) (m^j)' 1 T(F), 
9ijk = Pi (Pj ~ <%) (p k - o~ik ~ Sjk) {HiUjUkT 1 T(F), 

and so on, where 5ij = 1 if i = j and otherwise. Set 

-= J i MF {x a ) HF (x b )---dF{x). 



So, [ 4 a ] = m(l a ) of (GOP and by f^J), and 



2 



[ij] = U>i-lMj-iM2 - ^2 W-iMj+i + W+i - ViHi 

where E™..-i m fh-im = YT fh-im is defined in Example 5.2. 
By 

- 2T 1 {F) = 2d =T(1 2 ) = T{F) {2 < 1, 2 > + < 1, 1 > + < l 2 >} , (5.15) 
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where 



< 1,2 >=J^PiPj [ij] N V;\ 



»<3 



<1,1>=^( K ) 2 [J/] M r 2 , 



<i 2 >=X)a 



Other terms are calculated similarly. For example, C2, T 2 (F) and 82(F) are given by $3.6}) , 
fPU , and g^j in torms o/T(l 2 ), T(l 3 ) and T(1 2 1 2 ). j4feo 6y (CO)) to dXTT]) 



T (l 3 ) = T(F)l (pj - <%) (p k - 5 ik - S jk ) {mHjVk) 1 [ 

ijk 



ijk] 



+ 3 X> ( Pj - <%) Ohm)- 1 [ 21 ] + Y^PiN 1 [i] }. ( 5 - 16 ) 



and 



t(i 2 i 2 ) = t(.fk y^Pi (pj - %) (j>k - Sjk - Sjk) (pi - Sji - Sji - Ski) (jiiVjUkm) 1 

^ ijkl 



^Pi (Pj ~ Sij) (Pk - Sik - Sjk) (niUj^k) 1 G 



ijk 



^Pi(Pj - Sij)(p,iPj) 1 H i:j + ^2,PiHi x m (l 2 l 2 ) L (5.17) 

A A A J 



where 



Gijk = 2[lj] [I] +4 [12 < l i 2 fc ] , 

fly = 4 [1,12 2 ] + [?] g]+2[12 i 12 i ], 



1°2^1 C 2^ • • • 



so i/iai 



. . . 



ah 



[12ilj2 fe ] = (i) 2 lii-2AjA k - i BjjA k 

jk 



for 



A? ~~ Mi+i — jp-j-i^ii Bij — Lti+j-i — jnj-ip,i — Hi-ifij. 
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By (33p, 

3 

[i/fc] = -ijkHi-ifij-iUk-i^ + ^ iWi-iVj-i - W2) 

3 3 
~ (A*J+fc+l ~~ Mi+l^fc ~~ Mfc+l^i) + Hi+j+k — ^ifJ-j+k + ^^i^j^k; 

[i j] = ~(02JA«i-2^j-lM3 + (*)2A*i-2 i i+2-A*iA'2) 

+ 2 U>i-i(M»+i ~ - 2 ^ (Mi+j - MiMj - Mi-lMj+l) ! 

[lil2|] = (j) 2 {(-3i/ii_i/ij_i + - mHj-2) + 2^ + i/ij_i} 

2 

[12*12,-] = (i) 2 {j)2^i-2^j-2l4 - 2 ^2 i U)2l^ifJ'j-2fJ-2 

+2ij (fii+j-2^2 - Mi-iMj'-lA*2 + Mi^i) • 
^4feo [*] /or 2 < i < 4 and //j(l 2 l 2 ) are given fry i f5. 6]) - 1X77]) . 

Example 5.8 Consider Example 5.7 with T(F) = fir- Then 

T (l 2 ) /T(F) = p [ 2 ] /v 1 + (p) 2 /v 2 [#] , 

T (l 3 ) /T{F) = p^ 1 [ 3 ] + 3(p) 2 /^ 2 [ 2 J] + (p) 3 /V 3 [£] , 

T (1 2 1 2 ) /T(F) = V ll- X ll r (1 2 1 2 ) + (p)2^ 2 H rr + (p)3^- 3 Crrr + (p) 4 /V~ 4 [^ " 



Irrl 



Example 5.9 Consider Example 5.8 with T(F) = Set j3 r = // r /i 2 r ^ 2 . Tnen 

T(l 2 ) /T(F) = -2p+(p) 2 (/3 4 -l), 

T (l 3 ) /T{F) = -6(p) 2 (/3 4 - 1) + (p) 3 (As - 3/3 4 + 2) , 

r (1 2 1 2 ) /T(F) = 12(p) 2 - 4(p) 3 (/?4 - 1 + 2/3f) + (p) 4 (ft - l) 2 • 

-Ti(F)/T(F) = Ci/T(F) = -p + (p) 2 (ft - 1) /2, 

C 2 /T(F) = (p) 2 (5/2 - ft) + (p) 3 {/3 6 /6 - ft - ft 2 + 5/6} + (p) 4 (/3 4 - I) 2 , 
T 2 (F)/T(F) =p+ (p) 2 (4 - 5ft/2) + (p) 3 (2ft - 9ft + 7 - 6ft 2 ) /6 

r 

+(p) 4 (ft - 1) 2 /8 = X)(p)<^ sa ?/> 
i=i 

4 

S 2 {F)/T(F) = (p) 2 (7/2 - 2ft) + 

For p = 2 t/iis gwes T(.F) = ^ 2 , 

C!=ih- 3/4 = -At4 + 3/4 (5.18) 

C 2 = -2/x 4 + 5/4 T 2 (F) = -5/i 4 + 10/i 2 ., 5 2 (F) = -4/x 4 + 7/4 (5.19) 

TVoie t/iai C\, C 2 agree with /^(2 2 ) of Sukhatme (1944, P a 9 e 368). 
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The UE of [i\ has the form 

l22=(j2a i22 (F)n-*)/l[(l-i/n). 

So, {ai = ai22(F)} are given by 



,i=0 / i=l 



a = T(F) = fi 2 , 

Ql = -6T(F) + Ti(F) = -/X4 - 3/4 

a 2 = 11T(F) - 6Ti(F) + T 2 (F) = /x 4 + 3//|. 



We now present a second method for finding an UE of Oi^f - This method avoids 
computing {Tj(F)}, but derives the UE of the vector 

that is, for all products of a given degree p, directly from their first few coefficients {Cj}. 
Suppose T(F) has dimension d = d p . Then 

d = AjT(F), 

where Aj is a d x d matrix of integers and Ao = Irf, the identity matrix. So, 

a(n)" 1 T(F) 

is the UE of T(F), where 



a{n) = Y, A 
But this is known to have the form 



00 

* in 





p-i 

T n (F)=A,/JJ(l-i/n), (5.21) 

i=l 

where 

r b/2] ^ 
a. = I E B * n_i I T ( F )' 

where Bj is a d x d matrix of integers with Bo = 1^. So, 
b/2] fp-i \ 

i=0 U=l J 

= {l-Di^e + Da^e 2 -...} 

x {I d - Aie + (- A 2 + Af) e 2 + (-A 3 + A X A 2 + A 2 Ai - A? ) e 3 + • • • } , 



20 



where D\(p) = (2)2/2 and I2 2 (p) = (p)z(p — l/3)/8. So, the UE (|5.2ip is given in terms of 
{Ai,i<p/2}: 

Bo = Id, 

Bx = -D l (p)I d -A 1 , 

B 2 = D 2 (p)Id + Di(p)Ai - A 2 + A 2 , 

B 3 = -£ 3 (p)Id - Ai(p)Ai - Di(p)(-A 2 + A 2 ) - A 3 + A X A 2 + A 2 Ax - A?, 
and so on. The method also applies to obtaining an UE for 

where n = fi(F). A third method (for p < 8) due to Fisher (1929) is given in Section 12 
of Stuart and Ord (1987). Their Tables 11 and 10, pages 554-555 may be used to verify 
Examples 5.8 to 5.11 after some labor. 

Example 5.10 Consider Example 5.7 with T(F) = (/X4,/Lt 2 )'. So, S5. 20\) holds with p = 4 
and d = [p/2] = 2. 

By 115. S|) ; ( 15. /or /X4, Ci = — 4^4 + 6/x 2 and C 2 = 6//4 — 15/U 2 , in agreement with 
/i(4) on Sukhatme (1944, page 368). So, by [5~ffl) . (51$ 

Al=(" 1 4 _ 6 3 )^ A 2=(_ 6 2 " 5 15 )- 

50, 

Bi = -6I 2 - Ai = f ~j ~g ) , B 2 = 11I 2 + 6A1 - A 2 + Ai = 

So, UEs of /j,4 and a\ are fJ,4 n (F) and H22n(F), where 

3 

fMn(F) = {^4 + (-2^4 " 6/i|) n- 1 + (3/i 4 + 9/4) n~ 2 } / JJ(1 - i/n), 

i=i 

and 

3 

/W* 1 ) = {l4 + (~A*4 - 3mI) n _1 + (u 4 + 3^1) n~ 2 } / JJ(1 - i/n). 

i=l 

Ta6/e 5.5 gives i/ie relative bias of S np (F) as estimated from two runs of 60,000 simulations 
for p < 2 and i 7 normal and exponential. For n = 100 and p = 2, the estimates are poor: 
see Example C.3. For p = 3 the bias is zero. 

[Table 5.2 about here.] 
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Example 5.11 Consider Example 5.7 with T(F) = (^5,^3^2)'- So, &5.20\) holds withp = 5 
and d = [p/2] = 2. 

By \5.12i) . \5.13\) for /i 5 , C\ = — 5/is + 10/x 3/ u 2 and C2 = IO/X5 — 50^H2, in agreement 
with /i(5) of Sukhatme (1944, V a 9 e 368). By I15.15}) - (5.17\ ), for fi2f^3, 



T(l 2 ) = 2fi 5 - 16^2, T (l 3 ) = -24 M5 + 72^2, T (l 2 l 2 ) = 96/^2, 
giving C\= \i§ — 8/i3^2 C2 = —4^5 + 24^3^2- So, 

/ -5 10 \ , . / 10 -50 
Al = ( 1 -8 J andA *={_ A 24 

So, 

Bi = -IOI2 - Ai = 

and 



-5 -10 
-1 -2 



B 2 = 35I 2 + 10Ai - A 2 + A 2 - i H) 20 



1 5 

That is, UEs of fi§ and are fi§ n (F), and fi32n(F), where 

4 

M5n(-F) = {/^5 + (-5/^5 - IOM3M2) + (10/^5 + 20/i 3 /i 2 ) n~ 2 } / JJ(1 - z'/n 

i=l 

and 

4 

M32n(-F) = {^3^2 + (-A*5 ~ 2^3/^2) "T 1 + (/i5 + 5/i3^2) n -2 } / JJ(1 - i/n). 



i=l 



Example 5.12 Suppose k = s\ = 1 and T(i ? ) = g(fi2)- Sei g r = g^(fi2), and (3 r 
fj, r fi 2 r ^ 2 - Then 

fix = Hf{x) = X - fi, fl 2x = H2f{x) =[!%. — fi2, fi2xy = h^2F{x, fj) = ~2fl x fl y 

by By (AM>, 

\2\ =T(1 2 ) =5^22(1, 1)+5V2(1 2 ), 

where 

M22(l, 1) = J f4x = J f4x dF ( x ) = M4 - M2> 



^2 (l 2 ) = j H2xx = ~2fi2 by {OP- 
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Similarly, by &A.9\) to \A.11\) and \A.1S\) . 

9^222(1, 1, 1) + 3 5 2 M22(1, 1 2 ) + 3V2 (I 3 ) , 

S 4 M2222(1, 1, 1, 1) + 65 3 M222 (1, 1, I 2 ) + 9 2 {4/i 22 (1, l 3 ) + 3/i 22 (l 2 , l 2 ) } 
+5V2 (I 4 ) , 

5 V 22 (1, l) 2 + g 3 {2 M22 (1, 1) M2 (l 2 ) + 4^222(06, a, 6)} + (aV) 
+5 2 {4/^22 (a,a6 2 ) +M2 (l 2 ) 2 + 2/i 22 (a6, a6)| at a = b = I 



T(l 3 ) 
T (l 4 ) 



T fl 2 l 2 ) 



T(1 2 1 3 ) = 
and by jXTgjJ 



^a 4 a; say, 

1=2 



where 



r(lW)=J>^ 



1=3 



^222(1, 1, 1) = y /iL = Me - 3/U 4 /i2 + 2// 3 , 

M22 (l, I 2 ) = y ^2x^2xx = -2 (/i4 - /U 2 ) > 
M2 (l 3 ) = J H2xxx = 0, 

M2222(l, 1, 1,1)= y /4a; = A*8 - 4/X6/"2 + 6^4M 2 - 3 /4, 
M222 (l, 1, I 2 ) = y Ax^2xx = -2 (/i 6 - 2/U4M2 + Ml) > 



M 22 (l,l d ) =M2(1 4 ) =0, 

,2 



M 22 (i 2 ,i 2 )= y 



M2xx = 4/^4) 
M22 (a,ab 2 ) a=b=1 = J H2xfJ-2xyy = 0, 

H222 (ab, a, b) a=b= i 
/i 22 (a6,a6) a=b= i = 



[i-2xy^2x^2y 



M2 (a 2 6 2 ; 



0=6=1 



l^2xxyy — 0, 



a 2 = 12//1, a 3 = -4 (// 4 /m 2 - ^ 2 + 2/i 2 ) , a 4 = (/x 4 - ^l) 2 , 
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and 

u2\ , o.. ^„2\ ,. A ,2 



^3 = 6^222 (a, a&, & ) + 3/^2 (a ) ^22 & ) + 6^222 (b, ab, ab) at a = b = 1 

= ^J J {2fJ-2xV2xyH2yy + H2y V-2yy V-2xx + 2 //2j/^2x{/} 

= 3^ J {8(fJ? x - H2) ii x n\ + 12 (pjj - /z 2 ) M«A*J} 

= 12{2/x| + 3(/x 2 M4-^)} 

= 12^(2^1 + 3/34-3}, 

^4 = J J {^2xx^2y + Q^2xyf4y + ^^2yy^2y^2x} 

= -2 J J {nl (nl - fi 2 ) 3 + 6/i x ^ - nz) (nl - fi 2 ) 2 + 3f4 {nl - [I2) 2 {f4 ~ ^2)} 

= -2 1^2 (a*6 - 3/i 4 /"2 + 2/x|) + 6^3 (/x 5 - 2// 3 /x 2 ) + 3 (fj,4 - ^) 2 } 

= -2/4 {/3 6 - 3/3 4 + 2 + 6/3 3 (fc - 2ft) + 3 (/3 4 - l) 2 } , 

^5 = J J 'l4xl4y = J {L4-V2) 2 J (f4-^) 3 

= (M4 - M2) (A*6 - 3/U 4 /i2 + 2/i|) 

= Ml (& - 1) (& - 3/3 4 + 2) , 

S3 = 5* J ' fc at {a = b = c=l,S = n} 

= J J J {fJ>2xxH2yy/J>2zz + Q^2xx^2y Z + ^>^2xy^2yz^2zx] 



-120/& 



-B 4 = at {a = b = c=l,S = ^2} 

= 3 J J J {fJ^x^yy^zz + ^AxAyz + ^2x^2y^2xy^2zz + &H2xH2yH2xzH-2yz} 

= 36 { ( Ji 4 - M2) ^2 + 4 ^3>2 } 

= 36^{/? 4 - l + 4/3f} , 

-B5 = ^ j j j {^2xx^2y + ^2xy^2x^2y} ^2z 

= -6(^2(^4-^2) +^1} (/^-/i!) 

= -6^(/3 4 -l + /3 3 2 }(/3 4 -l), 

b 6 = / t4xt4 y f4z = (v4- 14) 3 = /AW*- 1 ) 3 ■ 
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So, 

Ci = -5V2 + g 2 (a*4 - iA) A 

C 2 = 5 2 (5^2/2-^4)+5 3 (^6/6-^3-^2 + 5^/6)+5 4 ( / u 4 -/U2) 2 /8, 

C3 = 0W2 + 5 3 (-Me/2 + + 2/4 - 6/4) 

+# 4 (ms/24 - /U 6/ u 2 /3 - // 5 /x 3 - f4/2 + 5/i 4 /4/ 2 + 5/4/*2 - 41/4/24) 
+0 5 (M4 " /4) (2a*6 " W2 - 3/4 + 7/4) /24 + 5 6 ( M4 - nl) 3 /48, 
Tx(F) = < S 1 (F) = -<7 2 (/"4-^)/2 + 5 V2, 

T 2 (F) = / (/i4 - /4) 2 /8 + / (/ia/3 -nl- 3/x 4 // 2 /2 + 7/4/6) + g 2 (-5/^/2 + 4/4) + ^2, 



where 



T 3 (F) = J>*T 3i , 
i=i 

5 2 (F) = 5 4 (M4-^) 2 /8 + /(w/3-^- 3^/U2/2 + 7^/6) + 5 2 (-2/m + 7/4/2) 

6 

5 3 (F) = 



i=2 



S32 = -3/X4 + 9/4/2, 

-533 = 3/t 6 - 27/t 4 // 2 /2 - 7/4 + 13/4, 

5 34 = -Ms/4 + 4// 6 /x 2 /3 + 2/X5/X3 + 11/4/8 - 6/M/4 - 7m§M2 + 85/4/24, 
«S 35 = (M4 - A*i) (-4/U6 + 15^ 4 ^2 + 3/x§ - 11/4) /24, 
536 = -5 6 /48, 
T 3 i = /is, 

T 32 = -19/^/2 + 31/4/2, 

T33 = 4/x 6 - 18// 4 /x 2 + 33/4/2 - 10/4, 

T34 = -Ms/4 + 4// 6 // 2 /3 + 2/X5/X3 + 7/4/4 - 27/^/4/4 - 7/4// 2 + 47/4/12, 
T35 = (^4 - (-4^6 + 15^ 4 /U2 + 3/4 - 11/4) /24, 
T 36 = -# 6 /48. 

Example 5.13 Consider Example 5.12 with T(F) = // 2 . Then 

9 1 = (9)*/4~'> 

T(l 2 )//4 = ( g ) 2 (ft-l)-2g, 

T(l 3 )/^ = (9)3 (#b - 3ft + 2) - 6(5)2 (ft - 1) , 

T(l 4 )/^ = (9)4 (ft - 4ft + 6ft- 3)- 12(g) 3 (ft - 2ft + 1) + 12(g) 2 ft, 

T(l 2 l 2 )//4 = ( g ) 4 (ft -l) 2 -4(9)3 (ft -1 + 2/4) + 12(9)2, 
T(1 2 1 3 )/^ = 12(9)3 (2ft 2 + 3ft - 3) 

-2(9)4 {ft - 3ft + 2 + 6ft (ft - 2ft) + 3 (ft - l) 2 } 

+(9)5 (ft - 1) (ft - 3ft + 2) , 
T(l 2 l 2 l 2 )//4 = -120(9)3 + 36(g) 4 (ft -1 + 4/3 2 ) 

-6(9)5 (ft - 1 + ft 2 ) (ft - 1) + (q)e (ft - I) 3 • 
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So, U = Ti(F)/T(F) and si = Si{F)/T{F) are given by 



ti = 8 1 = -(q) 2 (Pi-l)/2 + q, 

t 2 = (g) 4 (Pa - l) 2 /8 + (q) 3 (As/3 - 3/3 4 /2 + 7/6) + (g) 2 (-5/? 4 /2 + 4) + q, 
S2 = (g) 4 (Pa ~ I) 2 /8d + (g) 3 (ft/3 - $ - 3/3 4 /2 + 7/6) + (q) 2 (-2ft + 7/2) , 

6 6 

*3 = y^(g)^3i, S3 = y~l(g)iS3i 

i=l i=2 

for 

t 3 l = 1, 

i 32 = (31 - 19ft) /2, 

t 33 = 4ft - 18ft - 10$ + 33/2, 

t u = {-3ft + 16ft + 24ft ft - 84$ + 21$ - 81ft + 47} /12, 

t 35 = S35 = (ft - 1) (-4ft + 15ft - 11 + 3$) /24, 

% = S36 = -(/3 4 -l) 3 /48, 

«32 = -3ft + 9/2, 

s 33 = 3ft-27ft/2 + 13-7$, 

s 34 = {-6ft + 32ft - 138ft + 33$ + 85} /24 - 6ft - 7$ + 2ftft. 

Example 5.14 Consider Example 5.13 withT(F) = fi 2 , so ET(F) = (l-n~ l )T(F). As a 
check q = 1 above gives T(l 2 ) = -2/x 2 , T(l 3 ) = T(l 4 ) = T(1 2 1 2 ) = T(1 2 1 3 ) = T(1 2 1 2 1 2 ) = 
0, so h =t 2 = te = 1, si = 1, s 2 = S3 = 0. 

[Table 5.3 about here.] 

1/2 

Example 5.15 Consider Example 5.13 with T(F) = fi 2 = o~(F) say. Putting q = 1/2 
gives t\ = S\ = (ft + 3)/8, so an estimate of a(F) of bias 0(n~ 2 ) is 

cr^jl + n- 1 (ft(£)+3) /8}, 

where /3 4 (F) = ft = [x^^ 2 . To reduce the bias further use 
S2 = (l6ft + 22ft + 164 - 15$) /128, 
s 3 = ^240ft + 432ft - 2503ft + 2817 - 165$ 

+4764$ + 315$ - 560ftft + 420ft$ - 1920ftft^ /1024. 

Ta&Ze 5.5 gives the relative bias of S np (F) estimated from simulations for p < 2 and F 
normal and exponential. For n = 100 and p = 2, the estimates are poor: see Example C.4- 

The usual estimator of a(F) is the sample standard deviation, s.d. = {nfx 2 (F)/(n — 
l)} 1/2 , with mean ail-tln^+O^- 2 )}, where t\ = h-1/2. So, bias {s.d.} /bias {a (F)} = 
Ai + 0(n- 1 ), where Ai = (ft - l)/(ft + 3). 



26 



For the normal, exponential and gamma (*y), fa = 3, 9 and 3+6 7 -1 , so Ai = 1/3, 
2/3 and (57 + 12)/(6j + 12) and the s.d. improves on o~(F), although both are first order 
estimates, that is, both have bias 0(n _1 ). 

To see how S n 2(F) improves on the s.d., note that bias {S n 2(F)}/ bias {s.d. } = 
\2r1~ 1 + 0(n~ 2 ), where A2 = S2/t*. For the normal, exponential and gamma (7), 

fa = 15, 265 and I2O7" 2 + 1307 -1 + 15, 



.so 



s 2 = 65/64,767/32 and N(j)/64, 

A 2 = 65/16 rj 4.06, 767/32 rj 24.1 and N(X) (2.5 + 6A" 1 ) -1 /64, 
where N(j) = 6907~ 2 + 7887 -1 + 65. 

— 1/2 

Example 5.16 Suppose k = s\ = 1, T(F) = jj,/a = = g(fi, 1^2) = ft say. Again set 

fa = fj ir fi 2 r ^ 2 . Then the partial derivatives of g are g\ = [i 2 ^ ' > 9u = 0, 52 = -/i/i 2 3 ' 2 /^ 
512 = -/i2 3/2 /2, 522 = 3/i^2 5/2 /4, 5122 = 3 / u 2 5/2 /4, 5222 = -15 /^2 /8, an d so on. Set 
U^F) = n, U 2 (F) = fi 2 . Then defining lly ■ ■ ■ (l 7 , 1 J , . . .) as in UlJ^ - (AA^ , 

{7 U (1, 1) = J Uf x = J nl = fJ-2, 



U12O-, 1) = J U lx U 2x = J VxfJ-2x = M3, 

^22(1,1) = J Ul = fH-f4. 



So, by (AM, 



Also 



T (l 2 ) = fa + (3 fa + 1) /4. 



^122(1, 1, 1) = J Vxt4x = M5 - 2/U 2 M3, 
#222(1, 1,1) = J f4x = M6 - 3/^2 + 2//|, 
?7i2 (l, l 2 ) = y HxfJ>2xx = -2/i3, 
C/ 2 l(l,l 2 ) = J, 



^2x^xx — 0, 



So, by (HI 



#22 (l, I 2 ) = y ^2x^2xx = -2 (/i4 
#1 (l 3 ) = y Mxxx, 

^ 2 (i 3 )=y, to =o. 



T (l 3 ) /3 = (3/3 5 - 2/3 3 ) /4 + /3 (-5/3 6 + ll/3 4 - 6) /8. 
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Similarly, at (1, 1, 1 2 ), U 22 \ = 0, 

t/122 = -2 ((JL 5 - A*3A*2) , ^222 = -2 ()U6 - 2/U4//2 + , 

Lfy (1, l 3 ) = U % (l A ) = 0, U 12 (l 2 , l 2 ) = 0, U 2 (l 2 , l 2 ) = 4/i 4 , 

so 62/ JOS)) . 

T (l 4 ) = 3 (-5/3 7 + 3/3 5 - 3/3 3 ) /2 + 3/3 (35/3 8 - 132/3 6 + 242/3 4 - 97) /16. 
at a = b = 1, 

U 12 (ab,ab) = Ux {a 2 b 2 ) = U 2 {a 2 b 2 ) = 0, U 22 (ab,ab) = J \i\ xy = 4/4 
Ui 22 (ab, a, b) = U 22 (a, ab 2 ) = U\ 2 (a, ab 2 ) = U 2 \ (a, ab 2 ) = 0, 
U 222 (ab,a,b) = (J-2xyP2xfJ>2y = — 2^§. 



So, fry {Zjgp 

T(l 2 l 2 ) = 451222/W3 (/M - Ml) + 52222 (/i4 - - 45122^3^2 

-4^222 { (M4 - Ml) M2 + 2/i 2 } + 12^22^2 

= -3 (5/3 4 - 43) £3/8 + 3/3 (35/3 2 + 90/3 4 + 320/3 2 - 77) /16. 

So, 

5i(F) = T\{F) = -/3 3 / 2 - /3 ( 3/3 4 + 1) /8, 

S 2 (F) = (48/3 5 - 15/3 4 /3 3 - 23/3 3 ) /64 + /3 (-80/3 6 + 446/3 4 - 327 + 105/3 2 + 960/3|) /128. 
Note that T(1 2 1 3 ), T(1 2 1 2 1 2 ) and 83(F) may be calculated similarly using (A.ty . 

Note 5.1 In the one sample example above [i is the mean of X ~ F. In many cases 
Xj = h(Yj) ; where h : R — )■ i? s is a given transformation and Yi,...,Y n ~ G on 
i?* is the original sample. So, n(F) = J xdF(x) = j h(y)dG(y). Equivalently, we may 
replace fi(F) = j xdF(x) by fJ,(F) = j h(x)dF(x), so that /x x = h(x) — /x. Similarly, 
if s = 1 replace fJ> r (F) = f(x — fi) r dF(x) by f(h(x) — fi) r dF(x) so that ( f5.^| ) holds with 
hi = h Xi = h(xi) — fj,. A similar remark holds for several samples. 



The next four examples apply this idea to return times and exceedances. 



Example 5.17 Take k = 1, h(x) = I(x < a) for some a in R s , and T(F) = fj, . Since 
[i = F(&), T(F) is the return period of the event {X < a}, where X ~ F. But the 
case T(F) = was dealt with in Example 5.3 in terms of \i r . In this instance jjL r = 
fj, r (Bi(l,p)), where p = F(&), so fi 2 = pq, where q = 1 — p, ^ 3 = pq(l — 2p) and (14 = 
pq(l — 3pg). So, by Examples 5.6, 5.7 and Note 4-3 an estimate of the return period p^ 1 of 
bias 0(n ) is S n 4\p\ = 5^4 [p] ifp>l or l~ x if p < I, where < I < p, 

3 
i=l 

and Si[p] = Si(F) is given by Si[p\ = p^ 1 — p~ 2 , S 2 [p] = —p~ l + p~ 3 , S 3 [p] = 2p~ 1 +p~ 2 — 
2p~ 3 — p~ 4 . 
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The same formula with p = 1 — F(a) and p = 1 — -F(a) gives an estimate of bias 0(n~ 4 ) 
for the return time of the event {X > a}. Similarly, for the event {x E A} with p = F(A) 
and p = F(A). Similarly, we can apply Example 5.4 to obtain estimates of bias 0(n~ p ) 
for any smooth function g(pi, . . . ,pk) given independent nipi ~ Bi(rii,pi), 1 < i < k. This 
problem can also be solved by the parametric method of Withers (1987). 



Example 5.18 Suppose k = 1, X ~ F on R? and T(F) = Er(X) | (X G i), where A C R* 
is a measurable set, F{A) > and r : R l — ^ R is a given function. Then T(F) = ^2 = 
Hl(F)/n 2 {F), where ^(F) = f /i i (x)dF(x), /^(x) = r(x)J(x E A) and /i 2 (x) = J(x G A). 
So, {Tj, 5j, 1 < i < 3} are given in Example 5.2 in terms of the moments of h5. 1\) in which 
Xj t now needs to be replaced by hj^x). Set 

p = F(A), q = l-p, h= f (r(x) - dF(x). 

J A 

So, n[2 J ] = m(Bi(l,p)) is given for 2 < j < 4 in Example 5.17 and 

n [1*2*] = i iq > + (-^y (- P y g . 

Using I\ = simplification yields 

Sn4(F) = fnp- 1 {1 - tfVVfa " 1) + ?V7(n " 1)2 + ^V 3 (2p - l)/(n - 1)3} . 

Unlike Example 5.17, one does not need to know a lower bound for p, since fi\ = if p = 0; 
so, if p = one interprets Sn^F) as an arbitrary constant. This shows, surprisingly that 
the bias reduction problem for T(F) = Hi/p can be treated as a parametric problem, the 
parameters being (fj,\,p). The more general problem ofT(F) = g(ni,p) does not reduce to 
a finite parameter problem as it involves {J A r l dF, i > 1}. 



Example 5.19 The conditional distribution of exceedances is 

F u (x) = P(X-u<x\X-u>0) = {F(x + u) - F(u)} / {1 - F(u)} (5.22) 

for x > 0. This is fJ-i/fJ-2 with A = {y : y > u} = (u, 00), B — {y:x + u>y>u} = (u,x + u) 
and r{y) = I(y G B). So, Example 5.18 applies with [i\ = F(x + u) — F(u), fi2 = 1 — F(u). 



Example 5.20 The mean conditional exceedance is 

ix (F u ) = j xdF u {x) = /Ui//i 2 

for 

Mi = J(x- u) + dF(x), n 2 = l- F(u), 



where 

x + = 

So, r(y) = (y — and Example 5.18 applies 



x, if x > 0, 
0, ifx<0. 
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The central moments of F u of (|5.22p are not covered by Example 5.18 and are probably 
best dealt with by writing them as functions of the noncentral moments and applying 
Example 5.1 with /i = {J(x — u) l + dF(x) 1 i > 0}. A more direct approach is given by the 
following example. 



Example 5.21 Suppose T{F) = S{F U ) for F u of (5JB3) - Set C y (F) = F u {y). Then 
C v ((1 _ e )F + e5 x ) = F u (y) + eC y F {x) + O (e 2 ) , 



T((l-e)F + e5 x ) = S{F u (-) + eC F (x) + 0(e 2 ))=S(F) + e / S Fu (y)C F (x)dy + O (e 2 ) , 



Higher derivatives can be calculated from t5.23\) . 

Now let us apply the previous note with s = 1, t = r, h(y) = a'y, where a lies in R r . Set 
H = EY. Then the joint central moment ^i— r = E(Y — ■ ■ ■ (Y — fi) r is the coefficient 
of ai ■ ■ ■ a r /r\ in /x r (a'Y), so the same relation is true of their derivatives. The same is also 
true of the cumulants. This device allows us to derive results for multivariate moments and 
cumulants from their univariate analogs. 

For example, from Example 5.6, for a univariate random variable, fJ, 2 (x) = (x — fi) 2 — fi 2 
and ii 2 {x\,x<i) = — 2(a;i — \x){x 2 — fi). So, for a bivariate random variable, /ii2(x) = (x — 
/x)l(x - /x) 2 - /Ui2 and [i 12 (x 1 ,x 2 ) = -2(xi - //)i(x 2 - /x) 2 . 

We illustrate this device further with the problems of estimating multivariate moments 
and the correlation of a bivariate distribution and its square. 

Example 5.22 Suppose k = 1, s = 2 and T(F) = [i\ 2 . From Example 5.6 and the previous 
remark, an UE of \i\ 2 is fJ,i 2 /(l — at F = F. 

Similarly, we have 

Example 5.23 Suppose k = 1, s = 3 and T(F) = /U123. An UE of //123 is Mi 2 3/{(1 — 



n" 1 )(l - 2n~ 2 )} atF = F. 

Example 5.24 Suppose k = 1, s = 2, and T{F) = /^{mh/^} -1 / 2 ; the correlation of a 
bivariate sample. So, (HHP of Appendix A holds with S(F) = (^12,^11,^22) and g(S) = 
Si{S 2 S 3 )- 1/2 ■ We shall apply fOp . Set Vij ... = //,,...!//„//„•••; 1 2 . So, T(F) = v X2 . 
Now S^l 2 ) = f 5 ixx = -2 m , 5 2 (1 2 ) = f 5 2xx = -2/x u and 5 3 (1 2 ) = f 5 3xx = -2 m . 
Also g x = (^n/U22)~ 1/2 , 92 = -vu/nn, 93 = -^12/^22- So, giSi{l 2 ) = T(F){-2 + 1 + 1) = 
0. Similarly, 5i x = (x - f*)i(x - /i) 2 - M12, so 5n(l, 1) = / 5 2 x = H\\ 22 - Mi2> anc ^ 
similarly 512(1,1) = M1112 - MnMi2 3 5i 3 (l, 1) = M1222 - M12M22, (1, 1) = Mini - Mii> 
533(1,1) = /i2222 - M225 anrf ^23(1,1) = M1122 - M11M22- 5o, an estimate of bias 0(n~ 2 ) 
is T(F) - T(l 2 )/(2n) or T(F) - T(l 2 )/(2n - 2) at F = F , where by / COD . T(l 2 ) = 

^12(3^1111 + 3^2222 + 2^H22)/4 — V\\\ 2 — ^1222- 



and 



where C F (x) = // 2 l I{u < x < u + y) — Mi/i 2 2 I(u < x). So 




(5.23) 
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Example 5.25 Suppose k = 1, s = 2 and T(F) = /xf 2 {/xn/U22} 1 = ^i 2 j the square of 
the correlation of a bivariate sample. Again iA . 1\) holds with S(F) = (^12,^11,^22) but 
now g(S) = Sf(S 2 S 3 )-\ so 9l = 2T(F)S^\ g 2 = -T(F)S^, g 3 = -T(F)S^, g l% = 
2T(F)S~ 2 , gi2 = -2T(F)( < Si,S 2 )- 1 ) 913 = -2T(F)(S 1 S 3 )- 1 , and g 23 = T(F)(5 2 5 3 )- 1 . 
Again giSi{\ 2 ) = T(F)(-4+2+2) = 0. So, an estimate of bias 0{n- 2 ) isT(F)-T(l 2 )/(2n) 
or T(F) - T(l 2 )/(2n - 2) at F = F , where by [Ajfy, T(l 2 ) = 2vf 2 (v lul + u 2222 + 2v ll22 - 
2^1112 — 2^222) • 



6 Estimating covariances of estimates 

In this section, we give an estimate of bias 0(n~ 3 ) for V n (i ? ), the covariance of T(F), 
where now T(-F) is a q x 1 vector with components {T a (F), 1 < a < q}. After Example 6.1, 
we estimate the covariance of more general estimates of T(F). 

From the formulas for {K^ b } on Withers (1982, pages 66 and 67), 

00 

Vf{F) = covar (T a (F),T?(F)) = ^n^iff (F), (6.1) 



where 



i=i 



Kf{F) = = £ A a J J T $ ( I ) T§ ( * ) dn a {xy) 

= ^A a T^(a,a), (6.2) 

= e^e/ 3 ^(: a y )n[ a z )d. a{xyz) ,2 



a a b\ T e(b 
1 w x y J \ z 



+T ?{l x) T ^(y l)}dKa(wx)dK b (yz)/2 
EA«E^ (a 2 , a) /2 

+ £ A Q A a I £ (a 2 6, 6) + r^(ab, ab) j /2, 



— fa/3 + f /3a, 



T^(a,a) = Tp ( ° ) Tp ( ) dF a (x) 



1; / F I 3; 

T^(a 2 ,a) = yTj^ « ) 2| ^ J dF a (x), (6.4) 
T^(a 2 b,6) = JJt^I I b y y^ b y yF a (x)dF b (y), (6.5) 



(6.3) 
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and 



T a P(ab,ab) = J J T| 



a b 



T p 



a b 
x y 



dF a (x)dF b (y). 



(6.6) 



Also, setting V a ^{F) = K^(F) and differentiating, we have 



and 



a/3 I a 



a T , a 
x r \ x 



2 r 

T^(a,a) + J2 J n 



a a \ p I a 
V x ) F \y 



dF a (y), 



p I a a 

x x 



X 



a a 

x x 



F 

Tia 
F 



x y 



a a \ p I a 



y 



dF a (y) 



a a \ p a a 



x y 



+ fin 



x y 

a a a 
x x y 



dF a (y) 



a a 
x y 



p ( a 
F 



dF a (y) 



so that 
Ci (v a ^F 



z 

[T afl (a 2 , a) + T a/3 (a 2 b, b) /2} + 2T a/3 (a6, ab) - T a/3 (a, a) 



- b=a 



So, n" l Kf (F) given by (|fL2|) estimates V£ P (F) with bias 0{n~ 2 ) and n^K^F) + 



a/3. 



_1 r^aP', 



n 



2 L a P(F) estimates (F) with bias 0(n" 3 ), where 



L ap (F) = Kf(F)-C 1 [V ap ,F' j 

= £( A «- A ')E TQ, V> a )/ 2 

+ £ A a A 6 | T aP (A b) + T a ^ab, ab) j /2 



6=a 



If = 1 this reduces to 



L a @(F) = T ap (a,a) - 3T ap (ab,ab)/2 



at a = b = 1, so that 



(n-l)- l T a P(a,a)-3n-' 2 T ap (ab,ab)/2 



(6.7) 
(6.8) 
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at {F = F, a = b = 1} estimates Vn(F) with bias 0(n 3 ), where at a = b = 1, 

T a P(a,a) = J T%{x)T^(x)dF(x), 



and 



T^(ab,ab) = J J ' T$(xy)TP(xy)dF(x)dF(y). 



One may prefer to use n 1 — n 2 instead of (n — 1) 1 in (|6.8p . Remarkably, unlike the case 
k > 1, the estimate (|6.8p does not depend on T al3 (a 2 , a) or T a/3 (a 2 b, b) at a = b = 1. 

We now show how to estimate 

W„,(F) = cow T (n) (F), (6.9) 

where 

oo 
i=0 

is g x 1 and T = T. Clearly, T (n) (F) estimates T(F). Now 

W n (F) = Y, n-'-iWn (Ti, T,-) , 

»J>0 

where 

W n (T i5 T,) = awar (t«(F), T,-(F) 

has (a, /3) element 

(T i5 t,) = w n (if, r/ ) = y n 12 (F) 

of {EJ with (T 1 ,T 2 ) = (2f ,2/). So, 



^(F)=X;n-^[F], 



z=i 



where 



i+j+k=i 



and 



So, 



K k (T\T 2 )=K\ 2 {F) of ([EI]). 
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of and 

Kf[F}=Kf(F) + A a/3 , 

where 

2 

A a/3 = E ir i (r Q ,if 

and 

K\ (T a ,T^ =Kf(F) 

of ([El at = if. 

So, n^KfiF) and n" 1 ^^ + n- 2 L a ^[F] estimate W^(F) with bias 0(n~ 2 ) and 
0(n~ 3 ), respectively, where 

L a/3 [F] = Kf[F] - d (V Q/3 ,F) = L a(3 (F) + A Q/3 . (6.10) 



Alternatively, for k = 1, the sum of (JfTSJ and n~ 2 A a P at F = F estimates Wn^(F) with 
bias 0(n" 3 ). Now for p > 2, T np of flO) has the form T (n) of dSSD with T x given by (IP) , 
so that 

< ( : ) - -a. {** ( s ) - w + ( £ ) / 2 < 

and so 

^i(TMf) = -^A 2 {r^(a 2 ,a)+T^(a 3 ,a)}/2, 

<V] = E( A «- A ')E ra "(A 6 )/ 2 -E A 'E ra/3 ( a26 ' 6 ),= a / 2 
+ E AaAfe | E T ^ ( fl2& ' 6 ) + Ta ^ ab > aft ) } a 

v*m = E(v 2 - A2 )E rQ/3 ( a2 '°) 

+ E A - Afe i E TQ/3 (° 26 ' h ) + TQ/3 ( ab ' aft ) 1 / 2 



E a2 {e ra/3 ( a26 > 6 ) + 2TaP ^ ah ) - T °"^ a ) } 



For k = 1, at a = 6 = 1, this gives 

2 



= - E (a 2 , a) + T a ? (a 2 b, b) } /2, 

Kf[F} = T a,3 (ab,ab)/2, (6.12) 
awar (T^ p {F),Ti p (F)) = n~ l T^{a, a) + n- 2 T a/3 (a6, o6)/2 + O (n~ 3 ) 
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which, remarkably, does not depend on T(a 2 ,a) or T(a 2 b,b) to this accuracy - whereas 
L a P [F] does. 

Example 6.1 Consider again Example 5.1, that is k = 1, T(_F) = g(/x), where now g may 
be a vector {g a }. By [ZTfy - [ZHW at a = b = 1 

Kf(F)=T<*P(a,a)=g?g^[ij], 
T a \ab,ab)=g%g{ l ^k\ fi[jl], 
T a ? [a 2 , a) =g% £ f,[ijk], 

(a\b)=g« k gf^ijHkl], 

andKf(F), L a ^(F), Kf[F], L al3 [F] are given by fOj) . (6l\ ), fUlQ) . / fOT]) . [613) . Note 
that L al3 depends only on the first and second moments of F, even though depends on 
the third moments! 

Example 6.2 Consider Example 6.1 with g{ii) = ac' fj, / f3' fi = N/D, say, - that is, Example 
5.2. Since q = 1 we drop suffixes a, (3. Define fj,[-] and 5{ as in 115.1]) and 115. 3]) . Then at 
a = b = 1 

K^F) =T(a,a) = D- 2 fi 2 [SS], 

T(ab,ab) = 2fi 2 [5(3} 2 + 2/x 2 [5% 2 [/3/3] , 

T(a 2 ,a) =-2D- 3 f i 3 [65P], 

T (a 2 b, b) = 2D~ 4 {2h 2 [5P} 2 + // 2 [5% 2 } , 

where /U 2 [<5/3] = 5if3jfi[ij] and ^[aj3^\ = aiPjjk^i'ijk]- I n particular, for g(fi) = ^i//x 2; at 
a = b = 1 setting •jij... = ■ ■ ■ )fi^ 1 /ij 1 ■ ■ ■ , we have 

K^F) = T(a, a) = (/xi//i 2 ) 2 (711 - 2 7 i 2 + 722) , (6.13) 
T(ab, ab) = 2 Oi//i 2 ) 2 (711722 - 47i 2 7 22 + 2 7 | 2 ) , 

T (a 2 , a) = -2 (^i//i 2 ) 2 ( 7 u 2 - 2 7i22 + 7222 ) , (6.14) 
T (a 2 b, b) =2 (//i/yu 2 ) (2 7 f 2 - 5 7 i 2722 + 3 7 f 2 + 7 u 722 ) . 

Note that \6.13\) is in agreement with equation (10.17) of Kendall and Stuart (1977). 

Example 6.3 Consider Example 6.1 with g(fj>) = N p , where N = a'fj,, that is, we consider 
Example 5.3. In the notation there, with a = b = 1 

K 1 (F) = T(a,a)=p 2 N 2p a{2) , 
T(ab,ab) = p 2 (p - l) 2 N 2p a 2 2) , 
T{a 2 ,a) = p 2 (p-l)N 2p a {3) , 
T(a 2 b,b) = ( P ) 3P N 2p a 2 2) . 

In particular, for s = 1 and g(fi) = [x p , with a = b = 1 

T(a, a) = p 2 ^- 2 ^, T(ab, ab) = p 2 (p - 1) V^VI, 
T(a 2 ,a) =p 2 (p- 1)^-^3, T{a\b) = (pW^Vl- 
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For example, var{Jl x } or (if Note 4-3 needs to be applied), var{fi 1 I > /)}, where I > 
is a known lower bound for \fi\, can be estimated by 

f n2 = {n- l)- 1 /!" 4 ^ - Qn- 2 fL-^l 

or by 

T n2 I(\il\>l) 

with bias 0(n~ 3 ), where (/x, /X2) is (^,^2) at F = F. Alternatively, replacing n~ 2 in T n 2 by 
(n — 1)~ 2 and setting s 2 = fa n/{n — 1), the UE of fj,2, we obtain 

T* 2 = n^/TV - 6n- 2 Jl- e s\ T^I > /) 

as estimates with bias (9(n~ 3 ). 



7 Estimating the covariance of an estimate of bias 

The emphasis of this paper has been to reduce bias, not estimate it. However, a number of 
papers have given methods for estimating the variance of an estimate of bias for the case 
k = 1. See, for example, Efron (1981) and Davison et al. (1986). These papers provide 
bootstrap and jackknife methods of an order of magnitude less efficient computationally 
than the Taylor series method (also called the delta method or the infinitesimal jackknife 
when p = 2) used here. 

Suppose then T(F) is a q x 1 functional. Note that T(F) has bias n~ 1 B(F)/2+0(n~ 2 ), 
where B(F) = |2| = J2^aT(a 2 ). Its estimate n" 1 B(F)/2 has covariance n" 2 V(F)/4 + 
0(n~ 3 ), where 

v-hf) = E A «/ B ?(") B ?(°) dF <.w 

= E^{M::> s (::j-^wv) 

+ V f f T a ( aaa ( aa,S ) + [ [ [ T a ( aaa \ T 13 ( aaa \\ 

^ J J \ xxy J \yy J J J J \ xxz J \ yyz J j 

and dF a (x), dF a (y), dF a (z) are implicit in the integrals. Finally, n~ 2 V(F)/4 estimates 
covar |n _1 B(F)/2| with bias 0(n" 3 ). 

The same is true if we replace B(F) by B rap (F). If desired, one could apply Section 6 
to reduce this bias to 0(n~ 4 ). 

Note 7.1 In equation (2.6) of Davison et al. (1986) and the following line a factor 1/2 
should be inserted. So, the usual bootstrap and the usual jackknife estimates of bias as well 
as our estimate n _1 B(F)/2, all have bias 0(n~ 2 ). 
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Appendix A 



Here, we note and illustrate the following chain rule for the partial derivatives of 

T(F) = g(S(F)), (A.l) 

where S(F) is q x 1 and g : R q — >■ -R. 

First, suppose /c = 1, that is, F is a single d.f. Given r > 1, let s(y) : R T ^ R q be an 
arbitrary function. Set <9j = d/dyi. Then 

T F (xi---x r ) = ai---a rff (s(y)), (A.2) 

evaluated with s(y) replaced by S(F), and d\ ■ ■ ■ d r s(y) replaced by S F (xi • • • x r ). So, 
setting 

T\... r = Tp (xi, . . . ,x r ) , 

Su—r = SiF (xi, . . . , X r ) , 

gij... = didj • • • g(s) 
with 9j = d/dsi at s = S(F), we have 

T\ = giSn, T\2 = gijSnSj2 + giSiu, (A. 3) 

3 

^123 = 9ijkSilSj2Sk3 + gij Sn2Sj3 + <7i<Su23, (-A-.4) 

6 

^1234 = 9ijklSnSj2Sk3Sl4 + ^ S'jl£'j2<S'fc34 
/ 4 3 \ 

I X/ ^1^7234 + ^ S , ii2-Sj34 J + ffiSil234, (-A-.5) 

where summation over repeated suffixes ... is implicit, and by the multivariate version 
of Faa de Bruno's chain rule given in Withers (1984), for r > 1, 

r m(n) 

17T1 ' ' ' Si k7Zk , (A.6) 

k=l n 

where ^ m(n ) sums over all m(n) = r!/ n[=i (i!™ 1 "^!) partitions (7Ti ■ ■ ■ 7Tfc) of 1 • • • r giving 
distinct terms with rij of the 7r's of length i, and ^ n sums over {n € iV r , ni = k, ^ in; = 
r}. For example, 

3 

E S'il2'S , j34 = Sn2Sj34 + S'ji3S'j24 + Su4<Sj23- 

The reader can derive T123 from T12 using equation (2.6) of Withers (1983) to appreciate 
the labor-saving this rule gives. 

By [4c] of Comtet (1974) the general term can be written in terms of the multivariate 
exponential Bell polynomials, {-B^S)^...^}: 

r 

T\... r = s ^g il ...i k B rk (S)i 1 ... ik . (A.7) 
k=\ 
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This is a much easier form to use than (|A.6|) as these polynomials are immediately derived 
from the univariate polynomials B rk (S) tabled on pages 307-308 of Comtet (1974). For 
example, the table gives 

-Bzu(S) = 54, 
5 42 (S) = 45!5 3 + 35 2 2 , 
£ 43 (S) = GS 2 S 2 , 
i?44(S) = S 1: 



so 



-841(8)1! — 5'j 1 i 2 34, 

4 3 
-B42(S)j 1 j 2 = 5'i 1 l5'j 22 34 + 5i 1 i 2 S'j 2 34, 
6 

-B43(S)j 1 j 2 i 3 = ^ S^lSj^Si^i, 

and (|A?7l) for r < 4 reduces to ([0]) - ([0]l . 

Now suppose F consists of k d.f.s: the only change is to replace (xi • • • x r ) by 
wherever it occurs. So, in the notation of ([33]) . (|X3]) - (|A3|) imply 



a± ■ ■ ■ a r 
xi x r 



T (a 2 ) = g ij S ij {a 1 a) + giSi{a 2 ) , (A.8) 

T (a 3 ) = gijkSi jk {a, a, a) + 3^,% (0, a 2 ) + (a 3 ) , (A.9) 
T ( fl4 ) = 9ijkiSijki(a, a, a, a) + 6g ijk S ijk (a, a, a 2 ) 

+ 5 y {45^ (a, a 3 ) + 3,% (a 2 , a 2 ) } + <^(a 4 ), (A.10) 
T (a 2 b 2 ) = 9ijkiSij(a,a)Ski(b,b) 

+gijk{Sij(a,a)S k (b 2 ) + Sij(b,b)S k (a 2 ) + iSi jk (ab, a, b)} 

+ gij {2Sij (a,ab 2 ) + 2Sij (b,a 2 b) + $ (a 2 ) S,- (& 2 ) + 25^(06,06)} 

+ gi Si (a 2 b 2 ) , (A.ll) 

where 

(a 7 , a J , . . .) = y S lF (£) 5 iF (£) • • • dF a (x), (A.12) 

a r ) = a a with / columns, (A. 13) 

x J x x 

Sa (aV, . . . , a K b L , ...)= J ... J S iF ■ • ■ ) *V (g£ • • • ) dF a (x)dF b (y) , 

(A.14) 

and so on. Similarly, from (IA.7j) at r = 5 we obtain 

5 

T(a 2 6 3 )=^ 5n ... lfc ^-^, (A.15) 

k=l 
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where 

A i = S t (a 2 b 3 ) , 

A ij = 2Sij (a, ab 3 ) + 35y {b, a 2 b 2 ) + S t (a 2 ) Sj (6 3 ) + 6,% (ab, ab 2 ) + 3,% (b 2 , a 2 b) , 
A^ k = Sij(a, a)S k (b 3 ) + 3S ljk (6, b, a 2 b) + 6S ijk (a, b, ab 2 ) + 6S ijk (a, 06, 6 2 ) 

+35jfc (b,b 2 ) Sj (a 2 ) + 6Sij k (b,ab,ab), 
A» kl = Si (a 2 ) b, b) + 6% w (a6, a, 6, b) + 3^(6 2 , 6)S ifc (a, a), 

^4 1 5 = Si 1 i 2 (a, a)Si 3 i 4 i 5 (b,b,b), 

and from (|A.7p at r = 6 we obtain 

6 

T(a 2 b 2 c 2 )=^9i 1 -, k B^, (A.16) 



where 



fc=i 



B i = Si (a 2 b 2 c 2 ) , 

B ij = B[ j + B%+B%, 

3 

B[ j = 2^Sij (a,ab 2 c 2 ) , 

3 3 

B { 2 j = ^2 Si (a 2 ) Sj (b 2 c 2 ) + 4 ^ (ab, abc 2 ) , 

3 

B# = 2 ^ Sij- (a 2 6, 6c 2 ) + 4Sij(abc, abc), 
B^ k = Bf* + B$ k + B ."' 



3 ) 



- B i' ?fc = X^ 5 ' i ^ a ' a ) 5 ' fc ( fo2<:;2 ) + 4 $Z 5 ^' fc i a ,b,abc 2 ) , 

6 6 3 

I?*"" = 2 S'jfc (a, ac 2 ) Sj (6 2 ) + 4 Sy^ (a, ab, 6c 2 ) + 8 Sij k (a, be, abc), 

3 

Bp = SHa 2 )Sj (6 2 ) (c 2 ) +2^5, (a 2 ) 5, fc (6c, 6c) + 8S ljk (ab, be, ca), 
B ijki = B i^ki + B ^ki^ 

6 

B[ jkl = 2 ^2 S ij ( a2& » b ) S ki(c, c) + 8S ijk i(abc, a, 6, c), 

3 

B 2 1 = ^2 {Sij(a,a)S k (6 2 ) Si (c 2 ) + 2Sij(a,a)S k i(bc,bc) 
+AS ijk (a, b, ab)Si (c 2 ) + 8S ijH (a, b, ac, be)}, 

3 

£u-i5 _ ^ {Si 1 (a 2 )Si 2 i 3 (b,b)S i4 i 5 (c,c) + Si 1 i 2h {ab,a,b)S ii i b (c,c)} , 
B 1 6 = Si 1 i 2 (a, a)Si 3 i 4 (b,b)Si 5 i 6 (c, c) , 
and ^ m is interpreted in the obvious manner by permuting a, 6, c. For example, 

3 

S^- (a, a6 2 c 2 ) = S'y (a, a6 2 c 2 ) + iSy (6, 6c 2 a 2 ) + (c, ca 2 6 2 ) . 



39 



Similarly, if we now allow T and g to be r-vectors with components {T a } and {g a }, then 
by (TQ1) . T a P(a,a) of (|6.2p is given by 

T^-(a, a,...) = g?g? ■ ■ ■ S ir ..(a, a,...) (A.17) 

and T Q/3 '(ab,ab) of (|6.6p satisfies 

2 

T^(a&, ab) = <$0&Si*(a, 0)^,(6, b) + ]T 9? 9%S i]k {ab, a, 6) + g? S^ab, ab), (A.18) 

a/8 

where 

S ijk (ab,a,b) = J J s iF f 1 ^ ( x ) ^ ( y ) dFa ^ dFb ^- 
Similarly, (J63D, HESJI yield 

(a 2 , a) = {sg-Sj ifc (o, a, a) + gfS ik (a 2 , a) } 5 f, (A.19) 

and 

(a 2 6, b) = Igf^Sijia, a)S kl (b, b) + g% [Si (a 2 ) %(&, b) + 2%(a&, a, 6)] 

+gfSu{a\b)jgf. (A.20) 

Similarly, 

r Q ^(a6, a, b) = {g«S ijkl (a, b, a, b) + sfS^ab, a, 6)} g^gf . 
We now consider the case, where S(F) is bivariate, that is q = 2. Since Sij(a I ,a J ) = 
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T(a 2 ) = {giiS ll +2g l 2S 12 + g 22 S22}(a,a) + {g 1 S 1 +g 2 S 2 }{a 2 ) , (A.21) 
T ( fl3 ) = {giiiSiii + 3gn 2 Sn 2 + 3gi 22 Si 22 + g 222 S 222 } (a,a,a) 
+3 {sil^n + 512 (-512 + S21) + g 22 S 22 } (a, a 2 ) 

+ {siSi + 5 2S 2 }(a 3 ), (A.22) 

( a ) = {film 'S'llll + 4511125*1112 + 6511225*1122 + 4<7i222 5*1222 

+5222252222} (a, a, a, a) + 6{5m5m + 51125*112 + 2gi 2 i5i 2 i 

+522l5221 + 25i225l22 + 52225222} («, «, O. 2 ) 

+4 {5ii5n + 512(5*12 + S21) + 522 5 2 2} (a, a 3 ) 
+3 {5115*11 + 2g 12 S 12 + g 22 S 22 } (a 2 , a 2 ) 

+ {5i5i + 525 2 }(a 4 ), (A.23) 
T (a 2 6 2 ) = {5iiii5ii5n + 2<7ni25ii5i2 + 5ii225n522 

+25l21l5l25n + 4g(i2125i25i2 + 2#i2225l2522 
+5221l5225n + 252212 5225l2 + 52222 522 522} (a, 6) 

+{5iii5n5i + 2<7i2i5i25i + 522i5225i 

+5ii25n5 2 + 2gi 22 5i2 52 + 5222 5 2 2 5 2 } { (a, a) (b 2 ) + (6, b) (a 2 ) } 

+4 {5iii5in + 35n 2 5il2 + 35i2 2 5i 2 2 + 52225 2 22} (ab, a, 6) 

+2 {511S11 + 512 (5i2 + 5 2 i) + 52 2 5 22 } { (a, ab 2 ) + (b, a 2 b) } 

+ {5n5i5i + 512 (5i5 2 + 5 2 5i) + 5225 2 5 2 } (a 2 ) (b 2 ) 

+2 {g n S u + 2g 12 S 12 + 522 5 22 } (ab, ab) + {51 5 a + g 2 S 2 } (a 2 b 2 ) . (A.24) 

The convention here is that 

(g^S^ H ) (a 1 , ■■■) = 57r 1 5 vr2 (a 1 , . . .) , 

(g-KiS^S^ H ) (a 7 , . . .) (V, . . .) = 5^5^ (a 7 , . . .) S n3 (b J , . . .) . 

Similarly, for q = 2, splitting the third term in (IA.15D . gi^A l ^ k , into the six components 
corresponding to A^ k , the first is 

gijkSijk = {gilkSllk + 2gi2kS\2k + 522fc522fc} 

at (a, a, ft 3 ) and similarly for the second and sixth components. Similarly, for the three 
components of the fourth term, the first being 

r 2 \ 

gi - lSi-l = < y^^gijjjSijjj + 35ill25jii2 + 5il225il22 > 
[i=l J 

at (a 2 , 6, 6, 6), and for the fifth term 

gil-i 6 Si v ..i s = (5ll-5n_ + 25i2-5i2- + 522-522-) 

x (5-lll5-lll + 35-1125-H2 + 35-1225-122 + 5-2225-222) 

at (a, a,b,b,b), where 5 7r _5 7r _5_ 7r '5_ 7r ' is interpreted as g nn 'S nn i. 
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Similarly, for q = 2, the term B\ in (|A.16|) has the component 

2 

^gijSij = ky^guSu + 85126*12 
%=x 

at (abc,abc). The sixth component is 

+ 2gi2-5i2- + 922-5'22-) Q?-ll-5'-ll- + 29-12-8-12- + 9-22-S-22-) 
X (ff-llfi'-u + 25_i 2 S'_i2 + 9-22S-22) 

at (a,a,b,b,c,c), where 3 7ri _5 7ri _5_ 7r2 _5_ 7r2 _g_ 7r3 5 , _ 7r 3 interpreted as 0^ 7^3 7^3, and 
so on. 
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The nonparametric analogs of the terms for £2 and equation (D.l) of Withers (1987) needed 
for T2 and T3 - apart from those given in (|3.3p - (|3.5|) are as follows. Summation over a, b, c 
is implicit, where they occur. These terms are listed both for the purpose of checking and 
for application to other problems. Note that T2 requires 



22 
10 



|3| and 



22 
20 



-2A^|2| a 



and that T3 requires 

= A^{T(a 4 )-r(aV)}, 
-2XlT(a 3 ), 



23 
10 

23 
20 

2 2 2 



222 
110 



1 

2 2 2 

2 

2 2 2 

2 2 

32 
20 



1 231 



2 2 2 
2 

-2\lT (a 2 a 2 ) 



-2\ 2 a \ b T (a 2 b 2 ) , 



2 2 2 
12 



2 2 2 
2 10 



2A 3 T (a 3 ) for 1 < i < 3, 



4A 3 T(a 2 ) 



32 
10 



\ 3 a {T(a A )-3T (a 2 a 2 )} 



-6131. 



Also, 



23 




2 


2 


2 




2 


2 


2 


30 







3 










4 
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since K a (x\x 2 ), being quadratic in F a , has functional derivatives higher than two equal to 
zero. To illustrate the proof, 



222 
210 



X a X b X c d x lIp ( ^ ° ) d y V F ( ° ) (Ik c (z 1 z 2 ) T F ( a a b 
where U(F) = k XiX2 = K a (xix 2 ) and V(F) = k ViV2 = n b (yiy 2 ). Note that 



Vf 



unless c = b and 



unless b = c = a. Also 



b c 

yi zi 



and 



Up ( a a ) = - A n (*i) A ^ > 

f ° ) = A, ( yi A y 2 ) - Y, A * (l/i) Fa (2/2) , 
^ ' yiy2 



where A 3/ (x) = (F a (x))y = I(y < x) — -F a (:r). Integrate first with respect to x = (xix 2 ): 
since columns in Tp I ) are interchangeable we may replace Y^ 1X2 ^ ^ Since 



Tp( a a a )dF a ( Xi ) = 
\ x 1 x 2 y 2 J 

for i = 1, 2, and 

d x U (2/1 < &l) ^ ( z i < ^2)} = ^ [x\ - y\) 5 (x 2 - zi) dx x dx 2 

with 5 the Dirac delta function, 



(B.l) 



d x U F [ a a )Tp a a 

Vl zi J \ xi x 2 y 2 



-2T F 



a a a 



yi zi 2/2 



So, 



222 
210 



= -2Xl f 4 dn a (z lZ2 )T F ( a a a )d y V F ( a ) 
1 J V Vx Vi z i ) V z 2 J 



Integrate with respect to y = (2/12/2): (IB.ip implies the contribution from the last two out 



of the three terms in Vp 



is zero. Also, 



A z (?/i A y 2 ) = I (z < y x ) I (z < y 2 ) - F a (y x A y 2 ) 
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so 



dy A 2 (yi A 2/2) = 5 (yi - z) 5 (y 2 - z) dyidy 2 - 5 (yi - y 2 ) dy 2 dF a (y 2 ) . 



So, 



a a a 



y\ 2/2 zi 



dyVF 



Z-2 



a a a 

Z2 Z 2 Zl 



a a a 



Vi V\ zi 



dF a ( Vl ) . 



Now integrate with respect to z = (ziz 2 ): by (IB.ip the second out of two terms from 
dK a (ziz 2 ) contributes zero. So, putting 



l iF - 



z 2 )T F 



we obtain 

222 
210 



-2A^ J dF a {z 1 Az 2 )S.T F 



a a a 
yi yi z 2 



a a a 

z 2 Z 2 Zl 



0. 



T F 



a a a 
yi yi zi 



dF a (yi 



-2A 



T F 



a a a 

z z z 



dF a (z)- / dF a (yi) L>= —2\ a T [a 



Appendix C 

Here, we show how to estimate N, the number of simulated samples needed to estimate the 
bias to within a given relative error e. 

Note that T np (F) has bias —n~ p T p (F) + 0(n _p_1 ) and that S np (F) has bias — (n — 
+ Oin-P- 1 ) = -n~ p S p (F) + C^n^ 1 ). Suppose we estimate the bias of Y = 
Snp(F) by Z = Y - T(F), where F = iV -1 J2f=i Y h Y o = S n P (Fj) and F j is the empirical 
distribution of the jth simulated sample. Then EZ = ES np {F) — T{F) is the true bias of Y 
and we can write Z = EZ + (v n /iV) 1 / 2 {7V(0, 1) + o p (l))} as N -)■ 00, where v n = varYi = 
Vrn^ + Oin- 2 ) as n -> 00, and V7 1 = F T (F) = ^ A a T(a, a) with T(a, a) = / T F { a x ) 2 dF a (x). 
So, if Sp = Sp(F) 7^ 0, the relative error in the estimate of bias, 

(bias estimate - bias)/bias « — (v n /N) 1 ^ 2 J\f(0, l)n p S p (F) 

~ -Vr(F) 1 / 2 S' p (F)- 1 n p - 1 / 2 iV- 1 / 2 AA(0, 1) 

is bounded by a given number e with probability greater than 0.975 + O p (n~ 1 ^ 2 ) if 

2V T (F) 1 / 2 S p (F)~ 1 n p - 1 / 2 N- 1 / 2 < e, 

that is, if 

N > N epn = e- 2 n 2 P"V P , 

where (f> p = 4Vt(F)S p (F)~ 2 . This implies that for e = 0.1 and n large, say n = 100, it is not 
practical to carry out enough simulations to give meaningful estimates of bias unless p = 1 . 
This is reflected by the poor estimates of bias in the tables for the case p = 2 obtained for 
n = 100 using N = 10, 000. 
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Consider the following one sample examples. Set /3 r = /U r // 2 • For F = Af(0,l), 
Ha = 3, (jlq = 15, fis = 105 and for F = exp(l), [1% = 1, A*3 = 2, ^4 = 9, ^5 = 44, //6 = 305, 
/* 8 = 14,833. 

Example C.l Consider T(F) = Then Vr = [14 — fj^. Si = H2, 4>i = ^{@a — !)• So, for 
a normal sample <f>\ = % and fi2 = ^(F) needs 



For an exponential sample 0i = 32, so one needs four times as many simulations. Since 
82(F) = 0, 4>2 is not defined. 

Example C.2 Consider T(F) = /^|- Then Vt = 4/i| (a*4 ~~ /■*!) an ^ by Example 5.8, Si = 
— H4 + H2, S2 = — 4/i4 + 7/^2 so f or a un tt normal, Vt = 8, S\ = —2, (pi = 8, S2 = —29, 
02 = 0.1522 so iVo.iin = 800n and N .i2n = 152n 3 and for exp(l), Vt = 14,048, 5i = 30, 
01 = 62.44, S 2 = 87, 02 = 7.424, so JV .u n = 6,244n and iVo.i2n = 74.24n 3 . 

Example C.3 Consider T(F) = ^4. T/ien Vt = /is — M4 ~~ 8/X5//3, and oy Example 5.6 
or 5.10, Si = 2(2/^4 — 3/if), S2 = 3(4/i4 — 7/if)j so /or a unit normal, Vt = 96, £1 = 6, 
01 = 32/3, S2 = 15, 02 = 128/75, so iVo.n n = 1067n and A/o.i2n = 171n 3 and /or exp(l), 
V T = 14,048, 5i = 30, 0i = 62.44, S 2 = 87, 2 = 7.424 ; so N .ii n = 6,244n and 
N .i2n = 74.24n 3 . 

1 12 

Example C.4 Consider T(F) = a = fi 2 ■ Then Vt = ^2(^4 ~~ l)/4, so by Example 
5.15, for a unit normal, V T = 1/2, S\ = 3/4, 0i = 32/9, S 2 = 1/32, 2 = 2048, so 
AT 011n = 356n and iV .i2n = 204,800n 3 and for exp(l), Vt = 2, S x = 3/2, 0i = 32/9, 
S 2 = 213/8 = 26.625, 2 = 0.01129, so N .im = 356n and N .i2n = 1.129n 3 . 



Appendix D 

Here, we list the non-zero derivatives /J, r -l2—p = Htf{xi ■ ■ ■ x p ) for 2 < p < r < 6. They are 
obtained from (I5.4p in terms of hi = fj, x ., where = x — fj,, the first derivative of [i: 




80, OOOn simulations for e = 0.01 
800n simulations for e = 0.1. 
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M2-i = h - H2, 

M2-12 = -2hih 2 , 

M3-1 = hf — (13 — 3/11/X2, 

M312 = -3 (/if - /U 2 ) /i 2 - 3/ii (/i 2 - M2) , 

M3123 = 12/li/l 2 /l 3 , 

//4-1 =hf-fl4- 4/li//3, 

M4-12 = l2hih 2 H2 - 4 (/if - // 3 ) /i 2 - 4/ii - // 3 ) , 

M4-123 = 12 (/if - jU 2 ) /l 2 /l3 + 12/ll - M2) ^3 + 12/li/l2 (/l 3 - M2) , 
M4-1234 = — 72/li/l 2 /l3/l4, 

M5-1 = hl- fi 5 - 5/11//4, 

^5.12 = 20/ii/i 2 ^3 - 5 (/if - ^4) ^2 - 5/ii (/i 2 - /U4) , 

M5-123 = -60/ii/i 2 /i 3/ t/ 2 + 20 (/if - /x 3 ) /i 2 /i 3 + 20/ii (/i^ - /x 3 ) /is + 2O/11/12 (/13 - A*s) , 
M5-1234 = -60 (/if - jU 2 ) h 2 h 3 hi - 60h\ [h\ - fi 2 ) /i 3 /i4 - 60/ii/i 2 (/i 3 - jU 2 ) /14 

-60/ii/i 2 /i 3 (/14 - /i 2 ) , 
M5-12345 = 480/ii h 2 h 3 h 4 h 5 , 
M6i = h\- hq- 6/11/X5, 

M612 = 30/ii/i 2 ^ 4 - 6 (/if - /i 5 ) /i 2 - 6/11 (/i 2 - // 5 ) , 

M6-123 = -120/ii/i 2 /i 3/ u 3 + 30 (/if - /X4) /i 2 /i 3 + 30/ii (/i 2 - /X4) /i 3 + 30/ii/i 2 (/i 3 - /X4) , 

M6-1234/120 = 3/li/l 2 /l 3 /l 4 /X 2 - (/if - /i 3 ) /l 2 /l 3 /l 4 - /ll (/l 2 - jU 3 ) /l 3 /l 4 
-/ll/l 2 (/l| - /X 3 ) /l 4 - /ll/l 2 /l 3 (/ll - /X 3 ) , 

/"6-12345/360 = (/if - fj, 2 ) h 2 h 3 h 4 h 5 + /i x (h\ - // 2 ) h 3 h^h 5 + /ii/i 2 (/i 3 - // 2 ) /i4^5 
+h 1 h 2 h 3 (hi - fj, 2 ) h 5 + hih 2 h 3 hi (h\ - /x 2 ) • 

Note that 

r 

^.i2...r = (-l) r - 1 (r-l)r!n^', 

and 

r-l 

/Vl2-r-l = (-l) r (^/2) ^ (/if - /U 2 ) /i 2 • • • /ir-1, 

where sums ° ver all r — 1 like terms. 
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Table 5.1. Relative bias of S np (F) 
for T(F) = fi^ 1 estimated from two runs of 5000 simulations. 







n = 10 




n = 100 








p=l 


p = 2 


p = 1 


p = 2 


Norm (1/2,1) 


Run 1 


0.0773 


-0.0242 


0.0089 


0.0013 




Run2 


0.0916 


-0.0092 


0.0087 


0.0011 


Norm (1, 1) 


Runl 


-0.0780 


-0.0105 


-0.0149 


-0.0094 




Run2 


-0.0660 


-0.0040 


-0.0141 


-0.0087 


Norm (2, 1) 


Runl 


0.0208 


-0.0048 


-0.0046 


-0.0070 




Run2 


0.0202 


-0.0056 


-0.0056 


-0.0078 


Exp (1) 


Runl 


0.1096 


0.0120 


0.0052 


-0.0045 




Run2 


0.1062 


0.0184 


0.0062 


-0.0035 



Table 5.2. Relative bias of S np (F) 
for T(F) = estimated from two runs of 60,000 simulations. 







n = 5 

p = 1 p = 2 


n = 10 
p = 1 p = 2 


n = 100 

p = 1 p = 2 


Norm (0, 1) 


Runl 
Run2 


-0.3584 -0.1988 
-0.3572 -0.1947 


-0.1934 -0.0543 
-0.1871 -0.0460 


-0.0174 0.0021 
-0.0206 0.0012 


Exp (1) 


Runl 
Run2 


-0.4957 -0.2861 
-0.4943 -0.2851 


-0.2831 -0.0754 
-0.2964 -0.0923 


-0.0380 -0.0063 
-0.0399 -0.0082 



Table 5.3. Relative bias of S np (F) for T(F) = a. 







n = 5 

p = 1 p = 2 


n = 10 
p = 1 p = 2 


n = 100 

p = 1 p = 2 


Normal (0, 1) 


Runl 
Run2 


-0.1578 -0.0265 
-0.1592 -0.0277 


-0.0764 -0.0082 
-0.0745 -0.0080 


0.0281 -0.0045 
0.0003 0.0031 


Exp (1) 


Runl 
Run2 


-0.2278 -0.1019 
-0.2331 -0.1084 


-0.1251 -0.0422 
-0.1206 -0.0422 


-0.0158 -0.0029 
-0.0176 -0.0004 


Number of 
simulations / run 




10,000 


30,000 


30,000 
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